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INTRODUCTION 

It  is  widely  recognized  that  there  is  at  present  no  single 
mathematical  optimization  technique  superior  to  all  other  tech- 
niques in  handling  every  type  of  problem.   Every  method  has  its 
own  merits  and  shortcomings.   Consequently,  one  may  be  suitable 
in  solving  some  types  of  problems  but  becomes  cumbersome  In  solv- 
ing others.   Remembering  that  we  are  dealing  with  the  optimiza- 
tion of  a  process  and  that  "optimizing  a  process"  is  itself  a 
process,  we  would  be  absent-minded  if  we  forgot  to  optimize  what 
we  are  doing.   The  problem  now  facing  us  is  to  choose  the  most 
adequate  technique  to  solve  a  specific  type  of  problem.   For 
some  problems,  the  best  method  may  be  to  use  several  techniques 
jointly,  as  is  illustrated  by  Lee  (32) .   In  order  to  do  so,  a 
comparative  study  of  all  available  techniques  is  desirable. 

Since  many  of  the  processes  encountered  in  practice, 
especially  in  the  chemical  industry,  are  so  complicated  that 
finding  the  optimal  design  and  operating  plans  for  them  chal- 
lenges the  ability  of  the  best  engineer,  a  plausible  approach  to 
optimizing  such  formidable  systems  as  complex  chemical  plants 
and  processes  is  to  break  them  down  into  manageable  subsystems 
which  can  be  optimized  individually  and  subsequently  reassemble 
the  optimized  subsystems.   However,  the  difficulty  associated 
with  such  an  approach  lies  in  taking  proper  account  of  the  inter- 
actions between  the  subsystems,  because  policies  which  are  opti- 
mal for  the  separate  units  may  be  disastrous  for  the  ensemble. 
The  multi-level  system  theory  describes  effectives  ways  of 


decomposing  these  large  systems  into  component  subsystems.   The 
maximum  principle  is  a  very  powerful  method  in  solving  problems 
of  a  stagewise  nature.   Thus  we  restrict  our  discussion  to 
these  two  techniques. 

In  the  first  part  of  this  work,  we  present  a  review  of  the 
literature  on  the  multi-level  approach  of  process  optimization 
and  control  and  a  critical  examination  of  the  derivation  of  the 
multi-level  optimization  and  control  techniques. 

In  the  second  part,  we  discuss  briefly  the  discrete  maximum 
principle  and  extend  it  to  a  system  with  inequality  constraints 
of  a  completely  general  form.   We  also  propose  two  computational 
schemes  to  solve  the  so-called  two-point  boundary  value  problem. 
The  comparison  of  the  multi-level  approach  with  the  discrete 
maximum  principle  is  also  included  in  this  part. 

In  the  final  part,  we  develop  the  system  model  and  formu- 
late equations  for  reverse  osmosis  water  purification  for  the 
purpose  of  optimizing  the  process  by  means  of  the  multi-level 
approach  and/or  the  discrete  maximum  principle. 


PART   ONE 


A    STUDY   OF  THE   THEORY  OP  MULTI -LEVEL   SYSTEMS 


k 


CHAPTER  I.   INTRODUCTION 

Daily  activities  of  individuals,  of  gigantic  enterprises, 
and,  as  a  matter  of  fact,  the  activities  of  the  whole  economic 
system  of  a  country  have  been  based  on  a  simple  principle  of 
"optimality" .   Observation  of  the  similarity  between  the  problems 
encountered  in  engineering,  which  include  optimization  and  con- 
trol problems,  and  problems  which  arise  in  the  macro-economic 
theory  should  enable  us  to  derive  some  mathematical  techniques 
for  the  decomposition  of  large  optimization  problems. 

Large-scale  optimization  problems  have  been  a  constant 
source  of  difficulty  in  both  systems  engineering  and  operations 
research  since  their  inception.   Roughly  speaking,  an  optimiza- 
tion problem  is  considered  "large"  when  the  computational  re- 
quirement which  must  be  satisfied  in  order  to  find  the  optimal 
value  of  the  manipulated  variables  exceeds  the  capacity  of  cur- 
rent computing  machinery  or  when  the  quality  of  the  performance 
of  the  system  decays  significantly  in  the  time  required  to  com- 
pute a  new  control  solution  (31) •   By  adopting  basic  economic 
concepts,  we  should  be  able  to  develop  simultaneously  a  frame- 
work for  the  synthesis  of  "organization-like"  structures  and  at 
the  same  time  use  the  mathematical  interpretation  of  these 
"organizational  structures"  to  develop  efficient  computational 
algorithms  for  large-scale  optimization  problems  (31). 

In  a  multi-level  system,  the  overall  system  is  subdivided 
into  a  number  of  subsystems,  each  of  which  is  assigned  a  sub- 
objective  function  (goal)  and  a  performance  equation  (control). 


A  subsystem  may  represent  any  real  or  fictitious  entity  consist- 
ing of  a  finite  number  of  stages  (13). 

Another  level  of  control  is  assigned  the  object  of  coordi- 
nating several  of  the  subsystems  on  the  lower  level  and  these 
in  turn  are  coordinated  by  a  higher  level,  and  so  on.   The  re- 
sulting structure,  shown  in  Fig.  1.1,  is  triangular  in  form  with 
the  apex  ultimately  responsible  for  the  achievement  of  the  over- 
all system  object. 

The  conventional  problem  where  a  given  system  is  controlled 
in  a  manner  which  satisfies  some  pre-defined  objects  is  called 
a  single-level  optimization  problem.   By  "single  level"  we  mean 
that  in  general  no  managing  or  coordinating  controllers  are 
present.   A  characteristic  of  this  approach  is  that  although  the 
overall  system  may  consist  of  a  complex  of  interconnected  sub- 
systems, the  optimization  technique  cannot  take  cognizance  of 
this  fact.   As  a  result  the  solution  effort  usually  is  propor- 
tional to  the  square  or  cube  of  the  order  of  the  problem. 

A  multi-level  control  optimization,  problem  is  one  in  which 
the  structure  of  the  subsystem  is  acknowledged.   The  conven- 
tionally phrased  problem  is  subdivided  into  levels  of  organiza- 
tion, so  that  on  the  lowest  level  each  subsystem  can  be  optimized 
with  respect  to  a  subgoal.   The  subsystems  and  goals  are  co- 
ordinated at  a  third  level,  and  so  on  (5>)« 

The  major  advantage  offered  by  using  the  multi-level  ap- 
proach to  treat  optimization  problems  is  a  reduction  of  dimen- 
sionality which  is  especially  significant  for  large  systems. 
In  addition,  the  reliability  of  the  overall  system  is  not  limited 


3rd    level 


2nd    level 


1st  level 


System 


Fig, 


Multi-level     control     structure. 


S  =    Subsystem 
G  =    Goal 
C  =    Control 


The      arrows    represent    the     flow 
of     information    up     from    lower 
level      and     the     flow   of.   control 
signal      down     from     upper     level  . 


by  that  of  any  one  portion.   In  principle  the  subsystems  can  be 
arranged  so  that  failure  of  one  does  not  disrupt  the  perform- 
ance of  the  other,  since  they  are  operating  independently.   How- 
ever, the  overall  performance  will  be  affected.   A  conventional 
single-level  problem  probably  would  be  disrupted  by  failure  of 
a  subsystem  since  they  are  not  operating  independently.   Other 
additional  advantages  are  discussed  in  detail  in  references 

(13,  III-). 

The  price  we  must  pay  for  using  the  multi-level  technique 
is  the  cost  of  coordination.   Each  subproblem  is  solved  not 
once  but  many  times.   It  is  obvious  that  the  success  of  multi- 
level techniques  for  an  integrated  system  lies  in  the  decompo- 
sition of  the  system. 

CHAPTER  II.   LITERATURE  SURVEY 

The  concept  of  the  multi-level  approach  to  the  control  of 
interacting  systems  was  first  introduced  by  Mesarovic  and 
Eckman  (l),  Sprague  (2),  Sanders  (3),  and  Coviello  (14-).   The 
foundation  of  this  approach  is  to  distribute  the  effort  for  con- 
trolling a  system  among  several  sub-controllers  at  several  levels 

Later  Takahara  (3)  applied  the  multi-level  systems  theory 
to  linear  dynamic  optimization  problems  in  much  the  same  way  as 
Coviello  (I4.)  .   A  large-scale  control  system  which  performs  both 
the  optimization  function  and  control  function  is  decomposed 
into  small  subsystems  by  neglecting  the  interaction  between  sub- 
systems.  The  higher-level,  goal-seeking  units  compensate  for 
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the  neglected  interaction  through  successive  approximation  of 
the  intervention  parameters  as  suggested  by  Sprague  (2)  and 
Sanders  ( 3) • 

Lasdon  (6)  introduced  a  technique  for  the  discrete  process 
optimization  by  extending  the  concepts  developed  by  the  re- 
searchers mentioned  above.   He  treated  multi-level  problems  by 
introducing  a  pricing  mechanism,  an  ideal  long  used  by  econo- 
mists to  achieve  decentralization  in  economics  problems.   He 
designed  a  two-level  structure  which  allows  us  to  solve  a  class 
of  discrete  optimization  problems  by  iterated  solutions  of  sub- 
problems.   The  major  step  in  his  development  was  the  attachment 
of  prices  to  the  interacting  variables.   The  problem  was  solved 
by  iterating  the  prices.   These  prices  were,  in  essence,  the 
Lagrange  multipliers  of  the  integrated  problem. 

Pearson  and  Macko  (7)  extended  the  multi-level  systems 
theory  to  a  class  of  general  dynamic  optimization  problems  by 
drawing  on  some  of  the  ideas  presented  by  Lasdon  (6).   A  set  of 
intervention  parameters  to  decouple  the  subsystems  and  their 
goals  is  used  in  their  approach.   For  an  optimal  choice  of  the 
set  of  intervention  parameters,  which  is  determined  by  a 
higher-level  controlling  unit,  the  subgoals  of  each  first-level 
subcontroller  must  be  satisfied.   This  implies  the  satisfaction 
of  the  original  system  goal.   Pearson  (8)  also  considered  multi- 
level problems  from  the  variational  point  of  view. 

Brosilow  and  Lasdon  (9),  and  Lasdon  and  Schoeffer  (10) 
studied  some  multi-level  problems  using  the  classical  Lagrange 
method.   In  this  approach  the  second  level  uses  the  "price 
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adjusting  technique"  to  adjust  the  Lagrange  parameters  until  a 
set  of  parameters  is  found  such  that  the  subproblem  solutions 
solve  the  integrated  problem  (i.e.,  conventional  single-level 
problems) . 

There  are  several  papers  which  apply  the  multi-level  ap- 
proach to  control  system  design.   Lefkowitz  (11)  used  the  multi- 
level concept  to  break  down  a  large  overall  problem  into  sim- 
pler subproblems  as  follows:   (1)  the  process  was  decomposed 
into  subprocesses,  each  being  controlled  according  to  a  local 
suboptimal  performance  criterion,  and  (2)  each  subprocess  con- 
troller was  decomposed  into  a  hierarchy  of  control  functions 
which  distributed  the  load  and  responsibility  for  satisfying 
the  control  objective. 

Durbeck  and  Lasdon  (12)  presented  a  technique  for  objec- 
tively simplifying  complex  static  optimizing  control  models  by 
selecting  the  control  model  parameters  and  structure  to  maxi- 
mize performance.   They  showed  that  for  interconnected  systems 
of  high  dimensionality  the  resulting  parameter  search  may  have 
computational  difficulty.   A  two-level  decomposition  technique 
was  used  profitably  to  reduce  this  difficulty.   The  basic  par- 
ameter search  and  decomposition  techniques  were  also  used  in 
the  two-time  scale  control  approach  (11),  in  which  the  assumed 
structure  and  parameters  are  associated  directly  with  the  con- 
trol law  instead  of  the  system  models. 
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CHAPTER  III.   A  MULTI-LEVEL  FORMULATION  OF 
SIMPLE  FEEDBACK  PROCESSES 
(OR  STATIC  PROCESSES) 


1.   The  Conventional  Single-level  Optimization 
Problem  (or  the  Integrated  Problem) 


A  schematical  representation  of  the  simple  feedback  pro- 
cess is  shown  in  Fig.  1.2.   The  process  consists  of  N  functional 
subsystems  interconnected  in  series.   A  portion  of  the  output 
from  the  last  subsystem  is  fed  back  to  the  first  subsystem. 

The  nth  subsystem  produces  a  vector  of  finished  products 
yn'  (the  so-called  boundary  output  N  +  1  indicates  that  it 
goes  out  of  the  system) ,  and  a  vector  of  intermediate  or  state 

variables  xn,  which  serve  as  inputs  to  the  (n  +  l)th  subsystems. 

•  n  • 

It  receives  both  the  decision  variables  9   and  those  variables 

coming  from  (n  -  l)th  subsystem  xn   . 

The  steady-state  operation  of  the  process  is  described  by 

the  performance  equations. 

xn  =  Tn(en^  xn-l}  (l#1) 

yn,N+l  =  vn(0n,  x*'1)  (1.2) 

n  =  1,  2,    . . .,  N 
where  xn  is  an  sn-dimensional  vector  function,  yn>w+1  is  a 
^-dimensional  vector  function,  and  8n  is  an  con-dimensional 
vector  function. 

The  initial  feed  enters  the  system  at  a  rate  q,  whereas 
the  feedback  rate  is  r.  The  combination  of  the  feed  and  the 
recycle  stream  is  described  by  the  following  equation: 
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x°  =   M(xf,  x\  q,r)  (1.3) 


where  M  is  called  the  mixing  operator. 

When  the  flow  rates,  Q  and  r,  and  feed  stream  conditions, 
x^",  are  constant,  equation  (1.3)  can  be  rewritten  as 

x°  =  M(xN).  (1.1+) 

A  typical  optimization  problem  associated  with  such  a  pro- 
cess, neglecting  random  effects,  is  to  find  a  sequence  of  6n 
such  that  the  objective  function 

s  =  Z    fn(en>  yn'N+1)  (1.5) 

n=l 
is  maximized  or  minimized  subject  to  inequality  constraints 

Rn(6n,  yn'N+1,  x*1"1)  >  0  (1.6) 

n  =  1,  2,  .  .  . ,  N 
where  Rn  is  an  rn   dimensional  vector  of  function  9n,  yn*    , 
and  xn_1. 

The  problem  of  finding  an  optimal  decision  vector,  6n, 
which  satisfies  at  least  the  necessary  conditions  for  a  maximum 
or  minimum  will  be  called  the  integrated  problem  (5). 

We  assume  that  all  functions  defined  thus  far  are  at  least 
twice  dif f erentiable  in  all  arguments. 


■^Note  that  in  this  treatment  components  of  x*  must  be  con- 
sidered either  as  parameters  or,  if  they  are  free,  as  elements 
of  61. 
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2.   Multi-level  Approach 

View  each  subsystem  as  buying  and  selling  the  input  xn 
and  output  xn  from  and  to  other  subsystems  (6).   Associate  with 

each  xn  a  vector  of  prices  (real  numbers)  of  the  same  dimen- 

n 
sion,  p  ,  at  which  these  transactions  take  place.   Let  each 

subsystem  be  under  the  jurisdiction  of  a  manager  who  views  the 
inputs  x    as  being  independent  as  the  vector  0n  is.   This 
enables  us  to  separate  the  subsystems  by  cutting  the  relations 
between  them,  that  is,  the  performance  equations,  equation 
(1.1),  are  ignored.   In  doing  this,  we  break  up  the  overall 
problem  into  a  number  of  small  problems  ( subproblems) ,  each  of 
which  is  to  be  solved  by  a  real  or  fictitious  "first  level"  con- 
trol unit.   Thus  the  subproblem  is  also  called  the  first-level 
problem.   In  addition  we  synthesize  one  or  more  "second-level 
control  units  whose  function  is  to  coordinate  two  or  more  first- 
level  controllers. 

By  proceeding  in  this  way  we  hope  to  achieve  the  following 
economics.   If  the  process  is  a  real  one,  I.e.,  if  the  imagined 
organizational  structure  can  be  realized,  then  we  will  enjoy 
the  benefits  of  parallel  operation.   This  is  to  say  that  sev- 
eral parts  of  the  overall  problem  will  be  processed  simul- 
taneously. 

If  the  process  is  imaginary,  1. e . ,  if  it  is  simply  a  com- 
putational device,  then  we  have  traded  the  task  of  solving  a 
large  problem  for  that  of  solving  a  number  of  smaller  ones.   In 
either  case  this  procedure  may  lead  to  significant  computational 
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A.   The  Subproblem  (or  the  First-level  Problem).  Regarding 

pn  as  parameters,  the  nth  subproblem  is  for  subsystem,  n,  de- 
scribed by  the  relations 

yn,N+l  =  vn(Qn^  xn-l}  (l>2) 

Rn(6n,  yn'N+1,  xn_1)  >  0  ,  (1.6) 
in  which  a  set  of  9n  =  9n  and  xn  =  xn  that  extremizes  the  sub- 
objective  function 

sn  =  fn(0n^  yn,N+l}  +  (pn}T  Tn(0n^  xn-l} 

-  (pn"1)T  xn"!  (1.7)2 

is  found. 

The  subproblem  solutions  6   and  x   are,  of  course,  func- 
tions of  the  prices  pn.   Then  there  exist  values  of  pn  which 
are  designated  by  pn  at  which  the  subproblem  solutions  are  the 
optimal  solution  of  the  integrated  problem. 

To  prove  the  above  criterion,  let  us  define  the  Lagrangian 
for  the  integrated  problem  as 
N 
n=l 
and 


L  =  f;  (fn(en,  yn'N+1)  +  (pn)T  (Tn-  xn)  +  (un)T  Rn)} 

(1.8) 


p°x°  =  pNxN  (1.9) 


-'-Note  that  0n  denotes  the  solution  (or  optimal  value)  of 
the  original  integrated  problem,  and  6n  denotes  the  solution  of 
the  subproblems. 

(pn)  is  an  sn-dimensional  vector,  (pn)  is  the  transpose 
of  (pn),  and  (pn)T  T  (6n,  x11"1)  denotes  the  dot  product  of  sn- 
dimensional  vectors  (pn)  and  Tn. 


15 


where  the  vectors  pn,  n  =  1,  2,  ...,    N,  are  sn-dimensional  La- 
grangian  multipliers  for  the  equality  constraints,  equation 
(1.1),  and  the  vectors  un,  n.  =  1,  2,  ...,  N,  are  rn-dimensional 
Kuhn  and  Tucker  multipliers  (l6)  for  the  inequality  constraints, 
equation  (1.6) . 

By  using  equation  (1.9),  it  is  seen  that 

J-  ,    UsT     n    J-  /  n-lvT  n-1  , -,  -,nx 

Z-Cp/x=i(p    )   x    .  (1.10) 

n.=l  n.=l 

Substituting  equation  (1.10)  into  equation  (1.8),  the  La- 
grangian  equation  becomes 

L  =  E  {fn(6n,  yn)  +  (pn)T  Tn(0n,  x11"1) 
n=l 

-  (Pn"!)T  x*"1  +  (un)T  Rn}  (1.11) 

According  to  Kuhn  and  Tucker  (l6),  the  necessary  conditions 

for  an  extremum  are  that  there  exist  un  =  un  and  pn  =  pn  such 

that 

?L  dfn  m  <?Tn  rn   <?Rn 

=  +    (pii)l   +    (..n)1   =   0  (1.10) 

n        ^Qn  oan  3  Qn 


3Qn       dQn  d®  3Q 

+    (pn) (pn_1)T 


ax""1         ^xn_1  '      '       ^x""1 

-      T   ^ 
+    (un)T =   0  (1.11) 

un  >  0  (1.12) 

(un)T   Rn   =  0  ^  (1.13) 

-   on 


<?un 


=  R11  >  0  (l.llf) 
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—  =  xn  -  Tn(en,  x11"1)  =  o  (i.i5) 

for  all  n  =  1,  2,  . . . ,  N. 

All  quantities  in  equations  (1.12)  through  (l.l£)  are 
evaluated  at  the  optimum  point  6n.   But  the  conditions,  equa- 
tions (1.10)  through  (I.II4.),  are  equivalent  to  the  Kuhn-Tucker 
conditions  for  the  subobjective  function  Sn  subject  to  Rn  ^  0 
at  the  point  i 6n,  yn>    ,  xn  \   .   Thus  regarding  pn  as  prices, 
we  see  that  at  these  prices  the  subproblems  satisfy  the  neces- 
sary conditions  for  an  extremum  at  any  point  at  which  the  inte- 
grated problem  satisfies  the  same  necessary  conditions. 

B.   The  Second-level  Problem  and  an  Iterative  Scheme  ( 6) . 
The  fact  that  there  exist  prices  pn  that  decouple  an  integrated 
problem  is  utilized  to  derive  an  iterative  procedure  for  opti- 
mization. 

The  task  of  finding  the  optimal  parameters  p  ,  n  =  1,  2, 
...,  N,  is  delegated  to  second-level  units,  and  the  solution  of 
the  subproblems  for  a  given  set  of  P  is  the  responsibility  of 
the  first-level  units.    Note  that  since  all  the  subproblems  are 
independent,  the  first-level  units  need  not  communicate  with 
each  other;  i.e.,  the  constraints  given  by  equation  (1.1)  can 
be  ignored. 

With  the  subproblem  solution  6   and  x    substituted  into 
equation  (1.1),  we  can  get  the  supplies  from  the  nth  subsystem, 
i.e.,  T  (6  ,  x   ).   Then  from  the  difference  between  the  amount 


Recall  that  P  =  (p1:  p2 :  ...  :  pN)  . 
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xn  demanded,  which  is  determined  by  the  n  +  1th  subsystem,  and 
the  amount  Tn(6n,  x1-1   )  supplied,  which  is  determined  by  the 
nth  subsystem,  we  can  define  the  vector  of  excess  demand  for 
xn  as 

En(P)  -  xn  -  Tn(9n,  xn_1)  .  (1.16) 

It  is  evident  that  if  P  =  P,  En(P)  =0,  n  =  1,  2,  .  ..,  N. 
In  addition,  if  En(P)  =  0  for  all  n,  then  relations  given  by 
equations  (1.10)  through  (1.15)  are  satisfied,  which  implies 
that  P  =  P.   Thus  we  can  state  that  P  =  P  if  and  only  if 
En(P)  =  0  for  all  n. 

The  second-level  adjusts  the  parameters  P  by  a  price  ad- 
justment rule  suggested  by  Samuelson  (17),  i.e., 

d 

—  pn  =  En(P)  (1.17) 

dt 

With  the  excess  demands,  En(P),  formed  from  the  lower  levels  by 
equation  (1.16)  in  hand,  the  higher  level  can  apply  a  finite 
difference  approximation  to  the  price-adjustment  rule  given  by 
equation  (1.17),  which  will  yield  convergence  to  the  optimal 
prices  P  from  the  feasible  initial  guess  Pq . 

Now  the  operation  of  the  multi-level  scheme  can  be  de- 
scribed specifically.   Sequentially  it  proceeds  as  follows. 

1.  The  second-level  sends  to  the  first-level  units  an 
initial  set  of  parameters  p  ,  p  ,  ...,  p  . 

2.  Each  of  the  first-level  units  optimizes  its  sub- 
problem  using  these  parameters. 

3.  Inputs  and  outputs  of  the  first-level  units  are  trans- 
mitted back  to  the  second  level  which  forms  the  excess 
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demands  En(P) . 
l|_.  If  these  excess  demands  En(P)  are  nonzero,  the  second 
level  adjusts  the  parameter  P  by  the  price-adjustment 
rule  such  that  the  difference  will  be  reduced  and 
transmits  these  new  parameters  back  to  the  first- 
level  units. 
5-  The  process  is  repeated  until  the  excess  demands  are 

all  zero,  at  which  time  the  solution  is  optimal. 
C.   Convergence  of  the  Price- ad  jus tment  Rule .   It  has  been 
shown  (6,  13)  that  if  the  subobjective  functions  Sn  and  the 
constraints  Rn  for  all  n  are  concave  (for  maximization  problems) 
in  the  xn   and  the  6n  for  all  real  values  of  P  and  if  at  least 
one  of  these  functions  Sn  is  strictly  concave,  then  the  price- 
adjustment  rule,  equation  (1.17),  is  asymptotically  stable  in 

the  large  and  convergence  of  P  to  P  is  monotone  decreasing  in 
1 


The  stability  of  the  price-adjustment  rule  is  examined  by 
Lyapunov's  second  method   (29,  30).  A    Lyapunov's  function 


E   is  the  Euclidian  norm  and 


||  2 
E    is  defined  here  as 


1  I 

ETE.   ' 

Lyapunov's  second  method,  Theorem  II,  states  that  if  it 
is  possible  to  find  a  function  V(x)  which  has  the  following 
properties, 

V(x)  >   0,    for  x  f   xe  (equilibrium  point) 

V(x)  =0,    for  x  =  x~  (equilibrium  point) 
dV(x) 

and         ^  0,  except  for  the  possible  case  when  x  =  xp, 

dt 
dV(x) 

=  0,  then  the  system  is  asymptotically  stable. 

dt 
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dP 
chosen  for  —  =  E  is 
dt 


1 


V(P)  -  -i|E  ||2  (1.18) 


which  satisfies 


1..      ,2 


V(P)  =  -I  |E(P)  ||=0 
2  ' 

and  its  first-order  derivative  is 

dV      m  ^E  dP 

—  =  (E)T .  (1.19) 

dt        3?   dt 

dP 
By  substituting  —  =  E  into  equation  (1.19), 

dt 

dV      T  «2E 

—  =  (E)   —  E  .  (1.20) 
dt        3? 

dV 
Equation  (1.20)  shows  that  —  is  the  matrix  of  a  quadratic  form, 

dt 

By  Lyapunov's  second  method  it  is  seen  that  asymptotic 

stability  in  the  large  of  the  price-adjustment  rule,  equation 

dV 
(1.17),  requires  that  —  be  negative  definite  or,  equivalently, 

dt  2E 

that  in  equation  (1.20)  we  should  have  —  be  negative  definite 

3? 
for  all  P.   We  are  thus  led  to  consider  in  detail  the  elements 

2E  <?E 

of  —  .   In  Appendix  I  the  negative  definite  of  —  is  proved. 
^P  2P 

D.   Simple  Sequential  Process .   The  algorithm  derived  in 
the  previous  section  can  be  reduced  to  a  simple  sequential  pro- 
cess without  feedback. 

For  the  process  without  recycle,  the  ratio  of  feedback  r 
is  equal  to  zero,  and  equation  (1.3)  reduces  to 
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x°  =   xf  .  (1.21) 

It  is  necessary  to  remark  here  that  as  we  define  xn  as  an 
interstage  variable,     re  is  no  interstage  variable  entering 
the  first  subsystem  and  leaving  the  last  subsystem  (Nth  sub- 
system) .   x  and  x  must  be  equal  to  zero,  i.e.,  x  =  x  =  0. 
By  this  treatment  the  relationship  defined  by  equation  (1.9) 
will  be  automatically  satisfied.   It  is  worth  noting  that  in 
this  treatment  any  initial  conditions  must  be  treated  either  as 
parameters  in  the  first  subsystem  or,  if  they  are  free,  as  ele- 
ments of  0  .   Similarly,  fixed  right-end  conditions  must  be  in- 
corporated through  the  relations 

N,N+1    VN,_N    N-ls  (-     00, 

J    '  -   V  (9  ,  x    )  .  (1.22) 

3.   Extension 

Previous  sections  have  dealt  with  simple  sequential  pro- 
cesses.  In  this  section  we  consider  how  the  general  non- 
sequential systems  may  be  decentralized.   We  shall  see  that  the 
same  pricing  adjustment  scheme  will  suffice,  but  that  one  more 
parameter  is  now  required.   A  full  treatment  of  this  topic  is 
not  attempted.   We  merely  wish  to  indicate  that  an  extension  can 
be  made  and  to  demonstrate  some  of  its  major  features.   A  special 
and  yet  very  common  case  of  recycle  processes,  which  covers  many 
problems  considered  in  reference  (21),  is  treated  in  more  detail 
in  section  5>. 

A.   The  Integrated  Problem.   The  configuration  of  a  highly 
nonsequential  discrete  system  can  be  completely  described  by  the 
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equations  (18) 

xn  =  fn(yn,  0n)  (1.23) 

yn  =  wn(xn)  (1.21].) 

y0n  =   a0n  (1.25) 

xn,N+l  =  vn,N+l(yn^  Qn)                      (1.26) 

Dn(0n,  yn,  xn*N+1)  >  0  (1.27) 

n  =  1,    2,  .  .  .,  N 

where  xn£Rs   are  recirculated,  state  variables,   y  £R    is  the 
total  input  to  the  nth  unit  from  other  units,  and  x  = 
( x  ,  x  ,  ...,  x  ).   y   =  a   R    are  given  constant  vectors 
(that  is,  the  values  of  the  boundary  inputs  are  preassigned) . 
D  6.R    are  inequality  constraints,  0  £R    are  decision  varia- 
bles,  and  x  •»    6R       are  finished  products  (or  boundary 

outputs).   A  typical  unit  is  shown  in  Pig.  1.3. 

n 
The  optimization  problem  is  to  choose  a  set  of  6  ,  n  =  1, 

2,  .  ..,  N,  such  that  the  scalar  function  (the  objective  function) 

s  =  £T  Fn(en,  xn'N+1)  (1.28) 

n=l 

or,  by  combining  with  equation  (1.26), 

N 
S  =  ^7  Gn(0n,  y11)  (1.29) 

n=l 

attains  its  extremum  values. 

To  derive  the  optimization  algorithm  for  the  problem  we 

shall  first  assume  that  the  functions  fn(0n,  yn) ,  wn(x), 
vn,N+l(en^  7n)t    Dn(0n^  yn  xn,N+l}  and  Gn(yn  en}  arQ  contin_ 

uous  in  their  arguments  and  are  at  least  twice  dif f erentiable 
in  all  arguments.   Furthermore,  we  assume  that  a  set  of  optimal 
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decisions  denoted  by  0n,  n  =  1,  2,  .  ..,  N,  can  be  found. 

B.   The  Multi-level  System  Theory  Approach.   As  previously 
discussed,  the  original  overall  system  is  subdivided  into  a 
number  of  subsystems  each  of  which  is  assigned  an  optimal  sub- 
problem.   An  additional  level  of  problems  is  assigned  the  goal 
of  coordinating  several  of  the  subsystems  on  the  lower  level  and 
these  in  turn  are  coordinated  by  a  higher  level,  etc.   A  two- 
level  structure  is  treated  here. 

(i)   Formulation  of  the  first-level  problem 

For  a  subsystem  n,  described  by  the  equations 
xn,N+l  =  vn,N+l(yn^  en}  (1.26) 

and 

Dn(0n,  y11,  xn>N+1)  >   0  (1.27) 

we  are  to  find  a  set  of  0n  and  yn,  such  that  the  subobjective 
function 

Sn  =  Gn(6n,  yn)  +  (pn)T  fn(yn,  Gn)  -  (zn)T  yn     (1.30) 
attains  its  maximum  for  some  pn  and  zn  given  by  the  second 
level. 

(ii)  Formulation  of  the  second-level  coordination 
problem 

With  the  subproblem  solutions  6   and  y   in  hand, 
the  second-level  calculates  the  recirculated  state  variables  by 
equation  (1.23)  which  is 

-n  =  fn(jn  £n.)  (1.23) 

n  =  1 ,  2 ,  .  .  .  ,  N . 

With  calculated  'xrl,  a  new  set  of  the  parameter  zn  can  be 
adjusted  by  a  price-adjustment  rule 


2k 


dzn 

—  =  yn  -  wn(x^),    n  =  1,  2,  ...,  ...       (1.3D 

dt 

The  parameter  pn  can  be  computed  from  the  new  set  of  the  par- 


ameter zn  as 


N       Jw1 

Pn  =  r  (z^ 


(1.32) 
x=x 


i=l      ;>xn 
The  iterative  scheme  is  as  follows: 

1.  The  second-level  assumes  values  for  zn  and  pn  and 
sends  them  to  the  first-level  sub-problems. 

2.  Each  of  the  first-level  subproblems  optimizes  its 
subproblem  using  these  parameters. 

3.  The  first-level  solution,  6n  and  y  ,  is  send  to  the 
second  level,  which  forms  the  quantities  xn. 

l\..    If  equation  (1.31)  is  nonzero,  the  second  level  ad- 
justs the  parameter  zn  by  equation  (1.31)  and  computes 
the  new  set  of  parameters  pn.   These  new  parameters 
are  transmitted  back  to  the  first-level  units. 

$.    The  process  is  repeated  until  equation  (1.31)  is 
zero,  at  which  time  the  solution  is  optimal. 

The  extension  to  dynamic  problems  is  included  in  the  next 
chapter. 
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CHAPTER  IV.   A  MULTI -LEVEL  STRUCTURE  FOR 
LINEAR  DYNAMIC 'OPTIMIZATION  PROBLEMS 


Here  we  are  concerned  with  the  optimization  of  a  linear 
dynamic  system  with  respect  to  a  quadratic  objective  function. 
For  simplicity  we  assume  that  the  problem  is  stationary,  or 
time  invariant  in  the  specified  time  period. 

1.   The  Integrated  Problem 

For  a  system  described  by  a  linear  differential  equation 

(5) 

x(t)  =  Ax(t)  +  B6(t)  +  Cd(t),    0  <  t  1  T  (1.32) 

where  x(t)fRs  is  the  state  vector,  0(t)£Rr  is  the  decision 
vector,  and  d(t)^R   are  the  disturbances,  which  are  a  known 
function  of  time  t.   A,  B,  and  C  are  suitably  defined  constant 
matrices . 

The  boundary  condition  is  given  as 

x(0)  =  a  .  (1.33) 

The  problem  is  to  find  a  set  of  0(t),  0  £  t  <z   T,  such  that  the 
objective  function 

1    rT 

S  =  -  j   f(x  -  y)T  Q(x  -  y)  -  9Te}dt  (1.34) 

2  Jo 

attains  its  minimum. 


^By  a  quadratic  function  we  mean  a  homogeneous,  second- 
degree  expression  in  n  variables  of  the  form 

n    n 
F(x,y)  =  £f  '  "£T  a±Ax±-   yi)(xi  -  y-)  . 

f=l   j=l    1J    !      ^      J       J 
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y£Rs  is  a  reference  vector,  Q  is  constant  positive  sym- 
metric matrix,  and  xT  denotes  the  transpose  of  x. 

2.   Canonical  Equations 

To  obtain  the  solution  to  the  problem,  we  assume  that  the 
control  space  of  Rr  of  control  functions  9(t)  is  bounded  and 
twice  continuously  dif f erentiable  and  that  a  disturbance  space 
R^  of  disturbances  d  is  bounded  and  twice  continuously  differ- 
entiable.   We  also  assume  that  a  unique  set  of  6(t)  and  x(t), 
which  notes  S  minimum,  exists. 

By  means  of  the  calculus  of  variations,  we  define 

F  =  -  /(x  -  y)T  Q(x  -  y)  +  6Tej  +  zT(Ax  +  B6  +  Cd  -  x) 

2l  }  (1.35) 

Then  the  Euler-Lagrange  necessary  conditions  are 

—  p (__  F)  =  0  (1.36) 

d x        2 1  <?x 

and 

—  P  -  0  =  0  .  (1.37) 

Or  more  specifically,  we  have 

dz 
Q(x  -  y)  +  ATz  +  —  =  0  (1.38) 

dt 

and 

6  +  BTz  =  0  .  (1.39) 

At  the  point  where  t  =  T  the  transversality  conditions  re- 
quire that  for  all  ■  admissible  variations  dx,  d0,  and  dt  on  the 
surface  T  -  t  =  0,  we  must  have  (5) 
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dx  =  6x  +  xdt  ;      d6  =  69  .  ( 1 .  J4.O ) 

d  d 

(F  -  xT  —  P  -  6T  —  F)  ,  mdt 
dx  2  9    t_i 

+  (—  P)^=T  dx  +  (—  p)^=t  d9  =  0  .  (1.14-D 

dx  £9 

Since  69  and  6x  are  free  differentials,  the  above  condition 


becomes 

•T  ^       T1  d  ,  d        v  T1  •  /   ,  x 

F  -  x1  —  F  -  91  —  F+  ( —  F)1  x  =  9   at   t  =  T        (1.2+2) 

^x        <?9      <?x 

d 

—  F  =  9   at   t  =  T  (1.2+3) 
dx 

^ 

—  F  =  9   at   t  =  T  (I.I4I4-) 
P9 


j'U  -  y)T  Q(x  -  y)  -  9T9J  t=[p  =  9  (1.2+5) 


or 
1 
2 

z(T)  =  0  (1.2+6) 

6(T)  +  BT  z(T)  =  9  .  (1.2+7) 

Combining  equations  (I.38),  (1.39),  (1.2+6),  and  (1.2+7)  gives  the 
following  set  of  equations  called  the  canonical  equations. 

x  =  Ax  -  BBTz  +  Cd  (1.2+8) 

z  =  -Q(x  -  y)  -  ATz  (1.2+9) 

and 

z(T)  =  9  (1.50) 

where   9(t)  =  -BTz(t),    9  £  t  <  T  (l.£l) 

and    x(9)  =  a. 

Note  that  these  canonical  equations  can  also  be  derived 
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directly  from  the  maximum  principle  (18). 

In  order  to  apply  the  maximum  principle  let  us  introduce  a 
new  state  variable  xs+-^(t). 

1  A 
(t)  =  -  [      (  (x  -  y)T  Q(x  -  y)  +  eTe]  dt  (1.52) 

2  In     l  ) 


xs+l 

'0 


dxs+1(t)    1  „ 

x  ..  =  =  -  (x  -  y)1  Q(x  -  y)  +  e^  (1.53) 

3+1       dt       2 

then   xs+1(0)  =  0   and   xs+1(T)  =  S  .  (1.5k-) 

The  Hamiltonian  function  will  be 

H  =  zTx  +  zs+1  xs+1 

=  zT(Ax  +  B0  +  CD)  +  zs+1  -  Ux   -  y)TQ(x  -  y)  +  6Te|   (1.55) 

dz         2H 

—  =   z   =   -  -Aiz  -  z  +1  Q(x  -  y)  (1.56) 

dt        #x 

dzs+1  2H 

— 5±±  =  zs+1  =  - =  o.  (1.57) 

dt  ^xs+l 

The  boundary  conditions  are 

za+1(T)  =  1,      z(T)  =  0  .  (1.58) 

Substituting  the  boundary  conditions  into  equation  (1.57) 
we  have 

zs+1(t)  =  1,      0  <  t  <  T  .  (1.59) 

The  necessary  condition  for  H  to  be  an  extremum  with  re- 
spect to  6(t)  is 

dE 

—  =0  (1.60) 

29 

or 

1 
zTB  +  -  29  =  0  .  (1.60a) 

2 


29 


Thus  we  have 

0(t)  =  -BTz(t)  .  (1.61) 

Substituting  equation  (1.59)  into  equation  (l.j?6),  we  have 

z  =  -Q(x  -  y)  -  ATz  .  (1.62) 

The  performance  equation,  equation  (1.32),  can  be  derived  from 

the  Hamiltonian  function  as 

dx  ?E 
x  =  —  =  —  =  Ax  +  B6  +  Cd  . 
.  dt    dz 

Substituting  equation  (1.6l)  into  the  above  equation,  we 


have 


x  =  Ax  -  BBTz  +  Cd. 


(1.63) 


Thus  we  have  shown  that  the  canonical  functions,  equations  (l.£8), 
(l.6l),  (1.62),  and  (I.63)  which  obtained  from  the  maximum  prin- 
ciple, and  equations  (I.J4.8)  through  (l.£l)  which  obtained  from 
the  calculus  of  variations,  are  the  same. 

Now  to  apply  the  multi-level  multi-goal  structure,  let  us 
partition  the  performance  equation  and  the  canonical  functions 
into  N  £  s  subsystems. 
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f  z1 


I  . 
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zN 


Q: 


.    Q 
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5l   '    '    '   %t 
■    1  l\  T 


x1^1 
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(1.65) 


or 


e1(t) 


eN(t) 


N 


A 


N 


B: 


B1  ^  T 


.     .    B 


iW 


V 


z    ( t) 


xn   =  An   xn   +   Bn6n   +   Cndn 


n 


N 


n 

n      i 


n 


+     Z     (A?   x1    +   Bje1   +   C^d1) 
i=l 

zn  =   -Q^(xn   -  yn)    -    U£)T   zn 

-    E    {Q^x1    -    y1)    +    (A?)T    z1 


N 


T 


en  =  -(Bn)T  zn  -   Z    (B?r 

i=i    1 


xn(0)  -n 


(1.66) 


(1.67) 


(1.68) 


(1.69) 


i^n 
The  boundary  conditions  are 

a1 
zn(T)  =  0  . 
By  introducing  the  direct  and  indirect  intervention 


(1.70) 
(1.7D 
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variables,  sn  and  /n,  equations  (I.67)  and  (1.68)  are  simpli- 
fied to 

kn   =  A^xn  +  B^en  +  c£dn  +  E^  sn  (1.72) 

zn  =  -Q£Un  -  jn   +7n)  -  CA^>T  zn  (1.73) 

where 

n  — n.    -3-     ,    n.  i    n  — i    n   i,  ,    .  > 

En  s   =  Z   (Ai  x  +  Bi  9   +  Ci  d  )  (1.7W 

i^n 

Q£  /n  =  T    {q?  U1  -  y1)  +  (A^)1  z1  (1.75) 

i^n  ' 


3.   Multi-level  Multi-goal  Algorithm 

The  decomposition  of  the  original  integrated  problem  into 
the  two-level  multi-goal  problems  is  given  below. 

A.   The  First-level  Problem.   For  a  subsystem  n,  n  =  1,  2, 
...,    N,  described  by  a  linear  differential  equation 

xn  -  A*xn  +  B^6n  +  C*dn ■+  E^s11  (1.76) 

we  are  to  find  a  set  of  6n(t),  0  <  t  <  T,  such  that  the  sub- 
objective  function 
T 

sn  =  -  /    f  (xn  -  yn  +  rn)T  o£Un  -  yn  +  yn)  +  (en)T  en]<at 

2/0     l  '(1.77) 

attains  its  minimum  for  some  given  direct  and  indirect  interven- 
tion variables  sn(t)  and  yn{t)    with  the  initial  boundary 
condition 

xn(0)  =   an  . 
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The  solution  to  this  problem  provides  time  functions  xn(t) 
and  zn(t)  which  are  then  communicated  to  the  second  level. 

B.      Second-level  Problem.   Compute  a  new  set  of  sn(t) 
and  yn(t)  by  using 

E£sn  =  £    (aJx1  -  T     (B?)(b})T  J  +   C^d1}        (1.78) 

^n  n    *~~    /  n.  i     i.    /.n.T   it  /,  ^oN 

Qn^  =  Z-  JQi(x   -  y  )  +  (A±)   z  J  (1.79) 

i=£n 
with  the  xn(t)  and  zn(t)  computed  from  the  first-level  problems. 
The  operation  of  this  organization  is  as  follows  (5): 

1.  Initially  assume  yn   =  sn  =  0.   Then  each  subsystem  is 
assumed  to  be  independent  and  is  optimized  separately. 

2.  The  subsolutions  xn(t)  and  zn(t)  are  communicated  to 
the  second  level  from  which  the  proposed  intervention 
parameters,  yn   and  sn,  are  computed.   Here  xn(t)  and 
zn(t)  represent  the  proposed  time  history  of  resources 
and  price  levels. 

3-  Yn   and  sn  are  communicated  to  the  subsystem  and  each 
subproblem  optimized.   The  process  is  repeated  until 
it  converges  within  some  predefined  tolerance. 

Now  we  shall  prove  that  for  any  feasible  intervention  par- 
ameters, yn(t)  and  sn(t),  n  =  1,  2,    ...,    N,  the  subproblems 
have  unique  minima.   Furthermore  there  exists  an  optimal  inter- 
vention Yn{t)    and  sn(t)  such  that  the  subproblem  minima  co- 
incide with  the  integrated  problem  minimum. 

To  prove  this  let  us  notice  that  the  integrated  problem  and 
subproblems  are  of  the  same  class  for  which  existence  and  unique- 
ness are  guaranteed  by  a  positive  definite  Q,  and  hence  Qn. 


33 


Comparing  equations  (1.7&)  and  (1.77)  with  equations  (1.72) 
and  (1.73),  we  find  that  the  subproblem  satisfies  the  canonical 
form  of  the  integrated  problem  if  and  only  if  sn(t)  =    sn(t)  and 
yn(t)  =  /n(t)/  when.  6n(t)  =6n(t). 

The  principal  question  is  whether  yn(t)  and  sn(t)  converge 
to  yu(t)  and  sn(t)  respectively.  It  is  shown  that  there  exists 
a  T  > .0  such  that  the  multi-level  multi-goal  algorithm  produces 
a  convergent  sequence  of  intervention  functions  fyn,  sn T  with 
limits  |yn,  sn/  for  which  the  subproblems  solve  the  integrated 
problem  ( 5) • 

By  introducing  a  set  of  intervention  parameters,  which  are 
multipliers  or  prices  in  static  systems,  the  general  dynamic 
system  is  decomposed  into  a  collection  of  decoupled  small  sub- 
systems.  For  an  optimal  choice  of  the  set  of  intervention  par- 
ameters the  satisfaction  of  the  extremal  condition  of  each  sub- 
system implies  the  satisfaction  of  the  extremal  condition  of 
the  original  system.   Of  course,  the  conditions  which  guarantee 
convergence  of  the  intervention  parameters  are  rather  stringent 
when  one  considers  nonlinear  dynamic  systems. 


3i+ 


CHAPTER  V.   A  TWO-LEVEL  OPTIMIZATION  TECHNIQUE 
FOR  DISCRETE  NONSEQUENTIAL  PROCESSES 


Here  we  extend  the  multi-level  systems  theory  to  optimiz- 
ing a  class  of  general  highly  nonsequential  discrete  systems 
Lch  was  proposed  by  Fan,  et  al. ,  (18).   We  shall  make  use  of 
e  vector-matrix  notation  employed  by  Fan,  e_t  al. ,    in  the 
Continuous  Maximum  Principle  (18). 

1.   The  Integrated  Problem 

The  configuration  of  a  discrete  system  can  be  described  by 

the  equations 

xn  =  fn(yn  6n}  (l>8o) 

yn   =    r°na°n  +    T    rVnxV  (1.81) 

V  =1 

y0n   =   a0n  (1.82) 

yn,N+l  =  vn,N+l(yn^    en.)  {1.&3 

Ln(Gn,    yn,    yn>N+1)   >  0  (1.8U 

n   =    1,    2,     .  .  .  ,    N 
V>  =   1,    2,     .  .  .,    N 

where  6n£Ts  are  decision  variables,  xn£Rs  are  state  variables 

which  are  equal  to  the  recirculated  output  ynm,  m  =  1,  2,     ..., 
N,  yn£Rs  is  the  total  input  to  the  nth  unit,  which  is  a  linear 

combination  of  all  the  inputs,  and  yn^w+-i'iRs  is  the  boundary 
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1   On   s 
output  from  the  nth  unit,   y  cR      is  the  boundary  input  to  the 

nth  unit,  'a   cRs  is  given  constant  vector  (that  is,  the  values 

n   dn 
of  the  boundary  inputs  are  preassigned) ,  D  £R    are  inequality 

constraints,  and  P        are  diagonal  matrices  defined  as 


,Vn 


-  Vn 

0 
0 


oi,  r-) 

0 


0   0  .  .  .  .   0 
■vSn 


0  .  .  .  .   0 


0 


<A 


\>n 


,  n  =  1,  2,  .  .  .,  N 


1.85) 


yn 


where  the  diagonal  elements  -<r.  ,  k  =  1,  2,  .  .  . ,  s  are  non- 
negative  scalar  constants.   Equation  (l.8l)  means  that  the  state 
variable  of  a  unit  itself  is  the  output  variable  and  that  the 
kth  component  of  yn,  1  1L   k  <   s,    is  a  linear  sum  of  the  kth  com- 
ponent of  all  yun.   This  equation  is  the  boundary  condition  for 
the  unit  function  of  the  system,  equation  (1.80). 

The  optimization,  problem  is  to  choose  a  set  of  6n,  n  =  1, 
2,  . . . ,  N,  such  that  the  scalar  function  (the  objective  function) 


N 


nn/fin        n/N+1' 


Y;    Fn(Qn,    yIVi^x) 
n=l 


(1.86) 


attains  its  extreme  value.   Substituting  equation  (1.83)  into 
equation  (1.86),  the  objective  function  becomes 


A  stream  leaving  the  nth  unit  and  entering  the  V  th  unit 
will  be  identified  as  n\)th  stream. 
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s  =  2l   Gn(0n>  yn)  •  (1-87) 

nRL 

To  derive  the  optimization  algorithm  for  the  problem  we 
shall  first  assume  that  the  functions  fn(yn,  6n)  and  Gn(yn,  Gn) 
are  continuous  in  their  arguments  and  that  the  first  partial 
derivatives  exist  and  are  piecewise  continuous  in  the  argu- 
ments.  Furthermore,  we  assume  that  a  set  of  optimal  decisions 
denoted  by  6n  can  be  found. 

2.   The  Necessary  Conditions 

The  Lagrangian  for  the  overall  system  is 

L  =  £T  {an  +  (Pn)T  (fn  -  xn)  +  (un)T  Rn 

n=l  l 

+    (Xn)T    (     p0na0n   +    ^     rnVnxV    _   yn}j  (1>Q8) 

V>=1  ' 

where  pn  and  \n   are  the  Lagrange  multipliers  for  the  equality 
constraints,  equations  (1.80)  and  (1.8l),and  Un  are  the  Kuhn- 
Tucker  multipliers  (l6)  for  the  inequality  constraints  (I.8I4.). 
The  Un  are  constrained  to  be  non-negative  in 

Un  >  0  .  (1.89) 

We  can  rearrange  the  last  term  in  equation  (1.88)  into  a  form 
which  involves  only  the  inputs  and  outputs  of  a  single  sub- 
system (or  unit)  n,  as  follows: 

T    Un)T  (r0na0n  +  1l  rVn^  -  yn) 

n=l  tf=l 
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=  jr  (?vn)Tr0na0n  +  r  -^  uVr^V  -  ^  un)T  yn 

n=l  » =1    n=l  n=l 


xn 


-  r  un)    (r°na0n  -  yn)  +  E"  T   (^)  rnU 

n=l  n=l  V=l 

=    IT     {(^n)T    (P°na0n   -   yn)    +  T      (^)Tpn^xn]      .  ■  (1.90) 

n=l    L  V^l  I 

Therefore  equation  (1.88)  can  be  rewritten,  as 

l  =  Z   Un(Qn,   yn)  +  (pn)T  (fn(en,  yn)  -  xn) 

n=l  L 

+  (Un)T  Dn(0n,  jn,    yn'N+1)  +  Un)T  (P0na0n  -  yn) 


<-—  »n/nn   n   n    n  TTn   ,  1       ,N*  M  nn  n 

n=l 

where  ^n  is  tiie  part  of  the  Lagrangian  associated  with  the 
nth  subsystem. 

Necessary  conditions  for  a  maximum  of  the  objective  func- 
tion S  subject  to  constraints,  equations  (1.80),  (1.81).,  and 
(l.81|.),  are  that  the  Lagrangian  be  stationary  with  respect  to 
all  its  arguments  as  given  by  equation  (1.91).   These  conditions 
yield  a  set  of  vector  equations  which  must  be  satisfied.   They 
are  : 

—  =  p0na0n  +  V;  rVn^   -  yn  =  0  (1.92) 

d\n  v=i 
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—   =  —    fG"(en,  yn)  +  (pn)T  fn(6n,  yn) 

+  (Un)T  Dn(6n,  yn,  yn*N+1)}  =  0  (1.93) 

—  =  —  (Gn(en,  yn)  +  (pn)T  fn(en,  yn) 
2yn     <?yn  { 

+  (Un)T  Dn(0n,  yn,  yn'N+1)  -  (Xn)T  yn)  =  0      (1.9W 

—  =  fn(0n,  yn)  -  xn  =  0  (1.95) 

ah  n    ,  m   . 

=  -Pn  +  ^     U  )  Pni/  =  °  d-96) 

=  Dn(Gn,  yn  vn,N+l(en   n} }  >  0  (1>97) 

Un  >  0  (1.98) 

(Un)T  Dn(0n,  yn,  yn>N+1)  =   0  .  (1.99) 

n  =  1,  2,  .  .  . ,    N. 

It  is  worth  noting  that  equations  (1.97)  through  (1.99)  are 
the  Kuhn- Tucker  multiplier  conditions  (16)  for  the  inequality 
constraint,  equation  (I.8I4.). 


3.   Formulation  of  the  Multi-level  System  Approach 

A  complete  decoupling  of  the  subsystems  is  accomplished  by 
either  relaxing  the  subsystem  interconnecting  constraints, 
equation  (1.8l),  or  considering  the  state  variables  xn  as  var- 
iable parameters.   This  isolates  each  subsystem  from  other  sub- 
system^ and  creates  independent  problems.   In  the  second  decom- 
position method,  the  solutions  of  subproblems  always  satisfy  the 
integrated  system  equations.   It  is  called  the  feasible 
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decomposition  method.   In  the  first  method,  the  overall  system 
equations  are  not  satisfied  except  at  the  solution  of  the  inte- 
grated problem.   This  method  is  termed  the  nonfeasible  method. 
By  using  these  two  decomposition  methods,  the  original  inte- 
grated problem  can  be  decomposed  into  a  two-level  problem. 

A .   The  Algorithm  for  a.   Two-level  Problem  by  a_   Feasible 
Decomposition.    We  have  seen  in  equation  (1.91)  that  the  La- 
grangian  of  the  system  can  be  broken  down  to  a  collection  of  2n 
which  contains  only  variables  associated  with  a  single  unit. 
It  is  easy  to  see  that  the  necessary  conditions  for  the  La- 
grangian  to  be  stationary  with  respect  to  all  its  arguments  are 
identical  with  the  necessary  conditions  for  £      to  be  stationary 
with  respect  to  its  arguments.   This  gives  us  a  way  to  decompose 
the  whole  system  into  small  subsystems  or  units.   In  fact  the 
necessary  conditions  for  2n    to  be  stationary  with  respect  to 
its  arguments  are  also  the  necessary  conditions  for  the  subob- 
jective  function  of  the  subsystem,  n, 

sn  =  Gn(yn^  0n)  (1.100) 

to  have  an  extremum  in  6n,  yn,  subject  to  the  constraints 

xn  =  fn(yn^  en}  (1.80) 

n   /-i0n  On  ,  «c—  r-Ajn   tf  , ,  o-,  \ 

y  =  r      a    +  2_  P     x  (l.ol) 

DnOn,  yn,  yn'N+1)  >   0  (1.81^.) 


•J-In  this  method,  the  following  assumptions  are  made:  (a) 
the  vector  6n  has  at  least  as  many  components  as  the  vector  xn, 
and  (b)  the  Jacobian  matrix  (fn)  „  is  of  full  rank  for  all  6n 
and  x  . 
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In  this  formulation  components  of  xn  are  regarded  as  variable 
parameters.   This  suggests  a  two-level  structure  in  solving 
the  problem. 

(i)   The  first-level  problem 

For  a  subsystem,  n,  described  by  the  equations 
xn  =  fn(yn^  en)  (1.80) 

N 
yn  =pOnQOn  +  £~  pVn^  (1>8l) 

l)=l 

yn,x\T+l  =  vin,W(en)  yn}  (1>Q3) 

we  are  to  find  a  set  of  6n  with  Lagrange  multipliers  'pn  and 

\      such  that  the  subobjective  function 

Sn  =   c-n(en,  yn)  (1.100) 

attains  its  maximum  for  some  given  xn  and  Un.   Here  we  let  the 

inequality  constraints  be  considered  in  the  second  level,  that 

is,  the  Kuhn-Tucker  multipliers  Un  are  delegated  to  the  second 

level.   This  yields  unconstrained  and  therefore  simplified  sub- 

problems.   This  is  equivalent  to  a  modification  of  previous 

results. 

(ii)  The  second-level  coordination  problem 

The  second-level  adjusts  values  of  xn  and  Un  by 

using  the  Lagrangian  differential  gradient  method  (19),  that  is, 

dL 
(xn)i+1  =  (xn)i  +  k  (1.101) 

(Un)   _  k  if  Un  >  0 

(Un)i+1  =  j  (1.102) 

r"n  \  •  .p  rr'n  - 


(Un),  ,  if  Un  =  0 


kl 


and 

Dn(en,  yn,   yn'N+1)  >  o  (1.103) 

where 


Pn  +  r  rxv)irnV  -  o  (i.96) 


au 


n. 


=   Dn(0n,    yn,    yn>N+1)    >   0  (1.97) 


k  >  0  . 
The  "pn,  'Xn,  '6n,  and  y"11  are  the  solutions  of  the  first-level 
problem,  k  is  an  arbitrary  constant  which  is  to  be  so  chosen  to 
drive  equations  ( 1.101)  and  (1.102)  to  .convergence  as  quickly 
as  possible. 

A  computational  algorithm  for  this  two-level  structure 
then  proceeds  as  follows. 

1.  Assume  values  for  xn  and  Un  and  send  to  the  first- 
level  units. 

2.  Solve  the  subproblems  of  first-level  units.   This  yields 
6n(xn,  Un),  T^x11,  Un),  pn(xn,  Un),  and  yn(xn). 

3.  Send  the  first-level  solutions  6n,  A.n,  pv0 ,    and  yn  to 
the  second  level,  which  forms  the  quantities 

and  . 


9xn  3Vn 

2L        .        3L 

i}..  If  4   0  and  <0,    use  equations  (1.101)  and 

^xn  ,        ^Un 

(1.102)  to  generate  a  new  set  of  xn  and  Un. 

5-  Send  new  values  of  xn  and  Un  back  to  the  first  level 


and  iterate  from  steD  2  until  =  0  and  >  0. 

(l.lOij.) 

The  above  algorithm  is  guaranteed  to  converge  to  a 

global  maximum  of  the  original  objective  function 
S(61,  e2 ,  .  ..,  0N)  provided  S(61,  62,  .  ..,  6N)  and  Dn  are  con- 
cave functions  of  6  ,  6  ,  . . . ,  01  .   The  algorithm  converges  to 
at  least  a  local  maximum  of  S  provided  S  and  Dn  are  locally  con- 
cave.  The  details  are  given  in  Appendix  I  and  reference  (19). 

B.   The  Algorithm  for   Two-level  Problem  by  a_  Nonf easible 
Decomposition.   Now  we  let  the  conditions 

=  r0na0n  +  1C  PVnx  -  y11  =  0  (1.92) 

n  =  1,  2,  .  .  .  ,  N 

be  relaxed.   This  will  separate  the  subsystems  by  cutting  the 
links  between  them.   Then  assign  arbitrary  values  to  \nLR3    and 
U  £R   .   Equation  (1.96)  defines  the  Pn  in  terms  of  the  assumed 
values  of  \n,    equations  (1.93)-  through  (1.97)  defines  Gn  and  yn 


in  terms  of  P  ,  U  ,  and  A   and  gives  x   in  terms  of  9   and  y  . 
If  the  correct  values  for  the  \n   have  been  chosen,  then  the 
values  found  for  xn  and  jn   from  equations  (1.93)  through  (1.9?) 
will  also  satisfy  equation  (1.92).   If  the  Xn  and  Un  chosen  are 
not  correct,  -chen  equation  (1.92)  will  not  be  satisfied  and  we 
must  choose  a  different  set  of  multipliers  Xn ,  Un.   This  sug- 
gests a  two-level  structure  in  solving  the  integral  problems, 
(i)   The  first-level  problem 

For  a  subsystem  n  described  by  the  equations 
xn  =  fn(  n>  Gn}  (l#Qo) 


1+3 


yn,N+l  =  vn,N+l(Qn^  yn}  (1>Q3) 

find  a  set  of  0n\,  jn ,    xn  such  that  the  subobjective  function. 

s£  =   Gn(en,  yn)  +  (pn)T  xn  -  Un)T  yn  (1.105) 

attains  its  maximum  for  some  given  A.Xi  and  U  . 

n   — n 
It  is  evident  that  there  exists  ~k      =  X      such  that 

the  subproblem  solutions  solve  the  integrated  problems  (6). 

(ii)  The  second-level  coordination  problem 

n 
The  second-level  computes  a  new  set  of  \      by  using 

the  price-adjustment  method  (6,  l£)  . 

(Xn)  _.  -  Un)  .,  +  kEn  (1.106) 

l+l        i 

En  =  .=  p0nfl0n  +   ^  pvnx»    _   yn.  (1.107) 

^A.n  V7=l 

k   >  0 

n  =   1,    2,     .  .  .  ,    N 

where  k  is  an  arbitrary  constant  which  is  to  be  chosen  to 
drive  En  to  zero  as  quickly  as  possible. 

The  adjustment  of  all  or  some  of  the  Kuhn-Tucker  multi- 
pliers, Un\,  is  similar  to  the  feasible  method.   That  is, 

<9L 

/■(Un)  .  -  k  ,    if  Un  >  -0 

(UnJ1+1  -  J  3Un     '  (1.102) 

t  0  ,  if  Un  =  0 

and 

Dn(en,  yn,  yn'N+1)  >  o  .  (1.103) 

A  computational  algorithm  for  this  two-level  structure 

then  proceeds  as  follows. 

1.  Assume  values  for  \  ',   U  ',  n.  =  1,  2,     .  .  . ,    N,  and  send 

them  to  the  first-level  subsystems. 


kk 


2.  Solve  the  subproblems  of  the  first-level  subsystems. 
This  yields  solutions  pn,  %Ti{.vn,    Xn ,    Un)  , 

y  (p  ,  X  ,  U  ),  and  x  (0  ,  y  )  . 

3.  Send  the  first-level  solutions  to  the  second  level 

PL 
which  forms  the  Quantities  En  and  . 

£L 

.  .  If  En  r   0  and  *^  °>  use  eauations  (1.106)  and 

(1.102)  to  generate  a  new  set  of  X      and  U  . 

5.  Send  new  values  of  Xn   and  Un  back  to  the  first  level 

2L 
and  iterate  from  step  2,  until  En  =  0,  >   0,  for 

all  n. 

A  sufficient  condition  for  the  algorithm  to  converge  is 

that  each  subobjective  function  Sn  be  maximized  for  each  A.n  and 

Un.   Proof  of  the  convergence  is  given  in  Appendix  II  and 

reference  (6,  13) . 


Ij..   Discussion 

A.  Note  that  these  two  techniques  convert  an  optimization 
problem  of  high  dimensionality  with  inequality  constraints  into 
the  iterated  solution  of  a  number  of  smaller  unconstrained  sub- 
problems  and  simple  second-level  adjustment  procedures.   This 
represents  a  modification  of  the  results  in  section  3.   Since 
unconstrained  problems  are,  in  general,  much  easier  to  solve 
than  problems  with  inequality  constraints,  it  is  desirable  to 
delegate  the  Kuhn- Tucker  multipliers  Un  to  the  second  level 
since  this  yields  unconstrained  subproblems. 
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B.  In  the  nonfeasible  decomposition  method  we  can  inter- 

n 
pret  the  multipliers  \      as  the  prices  which  a  subsystem  must 

pay  in.  buying  its  feeds  jn   from  other  subsystems.   Similarly, 
pn  can  be  interpreted  as  the  prices  that  the  subsystem  charges 
for  its  products  which  go  to  other  subsystems.   The  subobjective 
function  St>  is  then  the  net  profit  a  subsystem  makes  from  all 
its  transactions  with  other  subsystems  and  with  the  outside  of 
the  system  boundary  with  this  interpretation. 

xn  becomes  the  supply  of  products  produced  by  the  nth  sub- 
system, and  yn  becomes  the  amount  which  the  nth  subsystem  de- 
mands.  The  function  En  is  then  the  excess  of  demand  over  sup- 
ply.  The  coordination  algorithm  represented  by  equation  (1.106) 
is  analogous  to  a  competitive  economy. 

C.  The  feasible  decomposition,  methods  have  been  used  by 
Mittem  and  Nemhauser  (23),  by  Aris,  Nemhauser,  and  Wilde  (2I4J, 
and  by  Nemhauser  and  Wilde  (25)  to  reduce  recycle  problems  to 
sequential  problems  which  can  then  be  treated  by  dynamic  pro- 
gramming.  It  can.  be  shown  that  feasible  methods  are  the  dual 
of  nonf  easible  methods  ( li|_)  .   Conceptually,  both  feasible  and 
nonf easible  methods  have  distinct  economic  interpretations, 
one  as  a  perfectly  competitive  economic  system,  the  other  as  a 
monopolistic  economic  scheme. 


PART  TWO 


THE  MAXIMUM  PRINCIPLE  AND  THE 
MULTI-LEVEL  SYSTEM  THEORY 


k-7 


CHAPTER  I.   INTRODUCTION 

Despite  the  recent  profusion  of  work  on  process  optimiza- 
tion, two  areas  still  pose  serious  difficulties.   The  first  of 
these  encompasses  problems  in  which  the  optimization  is  subject 
to  inequality  constraints  on  the  process  variables.   Such  con- 
straints are  often  necessary  to  make  the  optimization  meaning- 
ful, and  they  appear  in  a  variety  of  forms.   The  second  area 
deals  with  multidimensional  processes  of  complex  system.   Any 
optimization  is  here  made  difficult  by  the  implicit  nature  of 
the  process  model.   The  presence  of  inequality  constraints 
further  complicates  matters. 

In  this  part  we  shall  first  present  an  extended  version  of 
the  discrete  maximum  principle  for  the  optimization  of  simple 
staged  processes  subject  to  Inequality  constraints  of  a  com- 
pletely general  form  and  the  computational  schemes  to  solve  the 
so-called  two-point  boundary  value  problem.   Secondly,  we  shall 
discuss  the  interrelationship  between  the  discrete  maximum  prin- 
ciple and  the  two-level  structure  of  the  multi-level  system 
theory.   It  can  be  shown  that  the  adjoint  variables,  zn,  Intro- 
duced in  the  discrete  maximum  principle  will  have  the  same  func- 
tion as  the  Lagrange  multipliers  (called  the  prices),  p  ,  in- 
troduced in  the  multi-level  system  theory,  and  that  the  adjoint 

.    2E'n 

function,  i.e.,  zn   =  ,  is  a  necessary  condition  that  the 

^xn"! 
Lagrangian  be  stationary  with  respect  to  its  arguments,  xn   . 
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CHAPTER  II.   THE  DISCRETE  MAXIMUM  PRINCIPLE 


1.   The  Algorithm  for  the  Simple 
Feedback  Processes 


The  simple  feedback  process  described  previously  can  be 
represented  by  a  set  of  the  performance  equations 

xn  =   Tn(0n,  x11"1)  '  (2.1) 

yn,N+l  =  Vn(0n,  xn-!)  (2.2) 

n  =  1,  2,  .  ..,  N, 
the  mixing  equation 

x°  =  M(xN,  xf),   xf  =  given  (2.3) 

and  the  objective  function 

S  =  Z    fn(6n,  yn'N+1)  .  (2.k) 

n=l 

a  optimization  problems  associated  with  such  a  process, 

as  shown  in  the  previous  part,  can  be  solved  by  the  two-level 

system  theory.   It  can  also  be  shown  that  the  same  results  can 

be  obtained  from  the  discrete  maximum  principle. 

By  introducing  a  new  state  variable  xn  -,  ,  such  that 


with 


xs+i  =  x£i +  rn(e">  yn'N+1)  (2-5) 

n  =  1,  2,  .  .  .  ,  N, 


s+1 


Notice  that  in  order  to  obtain  clarity,  the  inequality 

constraints  Rn(6n,  x"n_  ,  yn>    )  >  0  are  temporarily  ignored. 
These  are  considered  at  the  end  of  this  section  where  it  is 
shown  that  none  of  the  main  results  are  altered  by  their  presence 
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we  can.  reduce  the  problem  to  the  standard  form  in  the  discrete 
maximum  principle. 

Equations  (2.1),  (2.2),  and  (2.5)  with  the  mixing  equation, 
equation  (2.3),  completely  specify  an  enlarged  process  with 
(  s  +  1)  state  variables  and  with  x  +n  as  its  objective  func- 
tion, i.e., 

S  =  T"1  c  xN  =   xN  (2  6) 

b  '  A,  cixi   xs+l  *  {d'0} 

i=l 

The  procedure  for  solving  such  an  optimization  problem  by 
the  discrete  maximum  principle  is  to  introduce  an  (s  +  1)  di- 
mensional adjoint  vector  zn  and  a  Hamiltonian  function  Hn 
satisfying  ( 21) . 

Hn  =  (zn)T  xn  (2.7) 


zn-l  =  . 


with  the  boundary  conditions 

s    n    2M,(xN 


2.8) 


Z7  -  \L     z°.   —J-^—   =   0  (2.9) 

X  A  —  I  J 


tr  z9  — i 

j=l   J    2x£ 

n  =  1,  2,  .  .  . ,  N, 


and 


i  =  1,  2, 


zs+l=1  (2'10> 


and  to  determine  the  optimal  sequence  of  the  decision  0n  from 
the  conditions 

=  0  (2.11) 

n.  =  1 ,  2 ,  .  .  .  ,  N , 


or,  if  6n  is  constrained,  the  optimal  decision  vector  6n  is  de- 
termined either  by  solving  equation  (2.11)  for  0n  when  6n  is 
interior  to  the  constrained  region  or  by  searching  the  boundary 
of  the  region  to  satisfy 

Hn  =  maximum  (2.12) 

n  =  1,  2,  ...,  N. 
It  may  be  noted  that  the  performance  equations,  equations 
(2.1)  and  (2.5),  can  be  written  in  terms  of  the  Hamiltonian 
function  as 

(2.13) 


(2.10) 


z^  =~^=  <+1  ,  n  =  1,  2,  ...,  N.       (2.14) 
s+1 


=  xn  . 

2zn 

In 

this 

pr 

oblem, 
s+1 

we  have 
=  1 

and 


Combining  these  two  equations  gives 


z£+1  =  1  ,      n  =  1,  2,  ...,  N.  (2.15) 


We  separate  the  (s  +  l)th  component  from  others,  that  is, 

T 


1 


n 


H"  =  (,«)   x"  +  z<<+1  xs+1 

=  (zn)T  xn  +  x":}  +  fn(en,  yn'N+1) 

=  (zn)T  Tn(9n,  x^1)  +  xn;l  +  fn(6n,  yn>N+1) 

3+1  (2.16) 


ind 


Note  that  _zn  and  xn  ape  truncated  vectors  with  s 
dimension. 
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z11'1  - =  (z  n)±  -= + .  (2.17) 

The  necessary  condition  fop  an  interior  maximum  is 

9En  „  *Tn(0n,  xn-!)    2fn 

=  (z^)1  = +  -  0  .  (2.18) 

The  optimal  solution  can  be  obtained  from  solving  simul- 
taneously equations  (2.1),  (2.17),  and  (2.l8)  with  the  boundary 
conditions  (2.9) . 


2.   Computational  Scheme 

As  shown  in.  the  preceding  section,  by  the  use  of  the  maxi- 
mum principle  the  optimization  problem  is  reduced  to  that  of 
solving  a  set  of  simultaneous  equations.   One  of  the  major  dif- 
ficulties in  such  work  is  the  solution  of  the  so-called  two- 
point  (or  split)  boundary  value  problem,  i.e.,  solving  a  set  of 
simultaneous  equations  with  mixed  boundary  conditions.   It 
represents,  in.  fact,  the  major  difficulty  and,  although  numer- 
ical analysts  have  given  considerable  attention  to  the  two- 
point  boundary  problem,  the  case  that  arises  in  the  study  of 
optimization  tends  to  be  particularly  difficult.   One  way  of 
solving  this  problem  is  the  so-called  steepest  ascent  iteration 
method1  (22) . 

It  is  recalled  that  according  to  the  maximum  principle, 


Sometimes  referred  to  as  the  "gradient  method  in  function 
space"  and  by  Merriam  (26)  as  the  "relaxation  method". 
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9n  '     jsen  to  maximize  the  Plamiltonian  Hn  as  a  func    i  of 

maximization  requires,  however,  knowledge  of  the  cor- 
rect zn.   Thus  there  are  two  ways  in  approaching  the  optimal 
point:   one  starts  from  guessing  the  decision  variable  9n,  and 
one  starts  from  guessing  the  adjoint  variables  zn  or  state 
variables  xn. 

Steepest  Ascent  of  the  Hamiltonian.   Suppose  that  the 
estimation  of  9n  is  corrected  in  such  a  direction  as  to  increase 
For  example,  at  every  time  the  following  rule  might  be  used 
to  proceed  from  the  ith  to  the  (i  +  l)th  approximation  (22) : 


e?,,,,  =  e*  +  k 


(i+D  ~  *(i)     ?Qn 


(2.19) 

en=e? 


i) 

where  k  is  a  suitable  positive  constant.   The  sequence  of  com- 
tations  for  one  iteration  would  then  be  as  follows: 

(a)  Given  the  ith  estimate  of  Gn  stored  in  the  computer 
(originally  a  guess),  computations  according  to  equa- 
tion (2.1)  are  carried  out  subject  to  the  specified 
initial  conditions,  and  the  resulting  values  of  xn 
are  stored. 

(b)  Given  the  stored  values  of  xn,  the  adjoint  functions, 
equation  (2.17),  are  calculated  in  reverse  stage  num- 
ber with  the  boundary  conditions,  equation  (2.9). 
These  equations  are  stable  if  the  transformation- 
equations  are  stable. 

(c)  At  each  step  of  numerical  integration,  as  zn  becomes 
available,  6n  is  up-dated  according  to  rule,  equa- 
tion (2.19) . 
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It  has  been  shown.  (22)  that  for  sufficiently  small  values 
of  k,  the  sequence  of  computations  for  this  iteration  would 
converge  to  the  optimal  point,  i.e.,  where 


<?Hn 


?Q 


n. 


=  0  . 


2.20) 


en=en 


B.   Adjusting  the.  Ad, joint  Variables .   Let  the  performance 
equation,  equation  (2.1),  be  relaxed,  which  means  cutting  of 
the  connection  between  units,  and  assume  values  of  zn,  n  =  1, 
2,  ...,  N.   Then,  from  equations  (2.17)  and  (2.18),  we  can  solve 
for  0n  and  x  ~  .   If  the  correct  values  for  zn  have  been  chosen, 
the  values  of  0n  and  xn   which  are  calculated  from  equations 
(2.17)  and  (2.l8)  will  satisfy  the  performance  equation,  equa- 
tion. (2.1).   If  the  values  of  zn  chosen  are  not  correct,  then 
the  equation.  (2.1)  will  not  be  satisfied  and  we  must  choose  a 
different  set  of  values  for  adjoint  variables,  z  ,  n  =  1,  2, 
.  .  .  ,  N. 

The  adjustment  of  the  adjoint  variables  depends  on  the 

amount  that  the  performance  equation  (2.1)  did  not  satisfy,  i.e., 

n        n       r   n.   mn,_n.   n-lx?         ,  , 

z(i+1)  ~  z(-\    +  k  jx   -  T  (6  ,  x    )j         (2.21) 

where  k  is  a  suitable  positive  constant.   The  computational 

sequence  for  this  Iterative  procedure  is  as  follows: 

(a)  Given  the  ith  estimate  of  zn,  n  =  1,  2,  ...,  N  stored 

in  the  computer  (originally  a  guess),  computations 

according  to  equations  (2.17)  and  (2.l8)  are  carried 

out  subject  to  the  boundary  conditions  (2.9),  and  the 

resulting  values  of  6   and  x   '  variables  are  stored. 


n. 


n-1 


b)  With  the  stored  values  of  0   and  x    equation  (2.21) 


Sk 


adjusts  a  new  value  of  zn. 
(c)  The  iteration  continues  until  equation  (2.1)  is 
satisfied. 


3.   Simple  Sequential  Process  With  Constraints 


or  a  process  with  constraints  on  both  decision  and  state 
variables,  as  given  in  the  rn  dimensional  vector  form 

Rn(Gn,  x11"1,  yn)  >  0  (2.22) 

the  necessary  conditions  for  an  optimum  can  be  obtained  by 
combining  the  Kuhn  and  Tucker  conditions  (lo)  with  the  algo- 
rithm of  the  discrete  maximum  principle  (21). 
The  Hamiltonian  is  defined  as 

Hn  =  (zn)T  x'n  +  (Un)T  Rn  .  (2.23) 

— n     — n 
The  necessary  conditions  for  a  saddle  point  in  9   and  u  , 

(un  2"  0),  that  is,  Hn  is  a  maximum  with  respect  to  Gn  and  a 

minimum  with  respect  to  un  ,  are  (l6) 

(2.1) 

(2.2^) 

num  (2.25) 

(2.26) 


xn  =   Tn(xn-i^ 

e 

'      3z" 

n-1 

B^1-1 

ZKn 

Hn   =  ma. 

—     U  .              OP 

3Qn 

un  >  0 

-MSiote  that  there  are  other  treatments  of  this  problem 
suggested  in  references  (18,  27,  28). 
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(un.)T  Rn  =  0  (2.27) 

and 

Rn ■>:   0  ,      n  '=  1/  2,  .  .  . ,  N  (2.28) 

with  the  boundary  conditions 


xf   = 

a 

zN   = 

1 

3=1     j  **? 

Tor  i  =  :s  +  1 

for  i  =  1,  2 ,     .  .  . ,    s 


(2.29) 
(2.30) 


Note  that  equations  (2.26)  through  (2.28)  are  the  so-called 
Kuhn-Tucker  conditions  (l6).  There  are  a  number  of  possible  al- 
gorithms for  solving  equations  (2.1)  and  ( 2 .  2l\.)    through  (2.30). 

One  computational  scheme  suggested  is  a  gradient  search 
on  the  surface  H  ,  which  simultaneously  ascends  in  6   (the  max- 
imizing variable)  and  descends  in  u   (the  minimizing  variable). 
The  technique  must,  of  course,  take  account  of  the  restrictions, 
u  —   0,    and  is  given  in  the  following  form  (19) : 

en  =  en  +  k  —  (2.31) 

«?Gn 
un  =  u n  +  k  un  (2.32) 

where  k  is  a  suitable  constant,  and 

r   0  ,     if  u?  =  0  and  R.  >   0 
•  n.    j    '  i  i  ,  , 

v±  =    j  (2.33) 

^  -R3?  ,   otherwise 
l 

J-    "~     J_  m  ^—9  9      *       *      9  * 

The  computation  proceeds  as  follows: 

n      n 

(a)  Choose  initial  values  for  6   and  u  . 

(b)  Find  the  corresponding.  xn  by  forward  solution  of 
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tie  performance  equation,  equation  (2.1),  with  the 

boundary  conditions,  equation  (2.29). 

(c)  Obtain  the  values  of  the  adjoint  variables  z  ,  by 

ba     -'d  solution  of  equation  ( 2 .  2L|_)  with  the 

boundary  conditions,  eauation  (2.30). 

Evaluate  and    .   If  these  values  are  nonzero, 

8n 
adjust  6n  and  un  by  equations  (2.31)  and  (2.32)  s 


return  to  step  (b). 


(e)  The  process  is  repeated  until  and  .u   are  zero, 

at  which  time  the  solution  is  optimal, 
be  above  algorithm  is  guaranteed  (15)  to  converge  to  a 
. obal  maximum  of  the  original  objective  function  S,  provided 
i rid  Hn  are  concave  functions  of  6  ,  n  =  1,  2,  ...,  N.   The 
algorithm  converges  to  at  least  a  local  maximum  of  S  provided 
S  and  Rn  are  locally  concave.. 

To  prove  the  above  assertions  it  is  necessary  to  prove 
gradient  of  the  Hamiltonian  without  inequality  con- 
straints is  the  same  as  the  gradient  of  the  objective  function 
S,  with  respect  to  6n.   Such  a  proof  can  be  found  in  references 
(30,  15) •   Once  we  have  the  desired  proof,  we  can  then  fall 
back  on  the  well  known  proofs  that  the  Lagrangian  differential 
gradie       Lod  converges  under  the  above  conditions  (29). 

Another  computational  scheme  would  be  to  start  by  assuming 
a  set  of  values  for  z  J  and  u  ,  n  =  1,  2,  . . . ,  N.   Then  compute 
forward  by  using  equations  (2.25)  and  (2.26)  together  with  the 
boundary  conditions,  equations  (2.29)  and  (2.30),  to  obtain  6n 
and  xn,  n  =1,  2,  ...,  N.   Using  0n,  xn  in  equation  (2.1),  we 
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can  check  if  the  assumed  values  are  correct.   if  not,  the 

assumed  values  are  improved  according  to  the  following  equations 

^Hn 
zn.  =  zn  +  k  (xn  _  ___)  (2.3I4.) 

n    n  ,  -.     •  n  /  0  0  ^  \ 

u   =  u   +  k  u  (2.35) 

where  k  is  a  suitable  positive  constant,  and 


0       if  un  =  0  and  Rn  :>  0 


un  =     .  (2.36 

•Rn     otherwise. 

Essentially,  this  algorithm  is  the  same  as  that  described 

in  the  previous  chapter.   And  it  has  been  proved  to  converge 

N   n      n 
asymptotically  to  the  correct  set  of  z   and  u  ',  if  such  sub- 
objective  function  is  at  least  locally  maximized  for  each  set 
of  zn  (6). 


CHAPTER  .III.   THE  INTERRELATIONSHIP  BETWEEN 
THE  DISCRETE  MAXIMUM  PRINCIPLE  AND 
THE  MULTI- LEVEL  SYSTEM  THEORY 


We  know  that  there  exists  a  close  relation  among  the  maxi- 
mum principle,  dynamic  programming,  and  the  calculus  of  varia- 
tions (l8).   In  this  chapter  we  shall  show  that  there  also 
exists  a  close  relation  between  the  maximum  principle  and  the 
two-level  structure  of  the  multi-level  system  theory. 

The  multi-level  system  theory  as  presented  in  Part  One  is 
based  on  the  decomposition  of  the  Lagrangian  of  an  integrated 
system.   The  necessary  condition  for  the  system  to  be  extremum 
is  that  the  Lagrangian  be  stationary  with  respect  to  all  its 
arguments,  which  include  the  Lagrange  multipliers.   By 
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decompos i       Lagrangian  of  the  original  integrated  probl<. 

le  system  is  decomposed  into  small  subsystems.   Then  a 
coordination  algorithm  manipulates  the  Lagrangian  multipliers 
of  the  subproblems  to  the  point  where  the  solutions  of  the  s 
problems  correspond  to  the  solution  of  the  original  integrated 
problem.   On  the  other  hand,  by  introducing  the  adjoint  vari- 
ables, the  discrete  maximum  principle  decomposes  the  overall 
extremum  condition  into  the  cascaded  extremum  conditions.   It 
has  been  shown  that  for  a  system  to  attain  its  local  extremum 
value,  it  is  necessary  to  choose  a  set  of  decision  vectors  such 
Hamiltonian  for  each  unit  is  stationary  or  extremum. 
srefore  if  we  can  show  that  the  Lagrange  multipliers  in  the 
Lti-level  system  theory  are  the  same  as  the  adjoint  variables 
in  the  maximum  principle,  and  that  the  adjoint  functions,  i.e., 


2xfl~1 


are  exactly  the  same  as  the  necessary  conditions  that  the  La- 
grangian be  stationary  with  respect  to  its  arguments,  x   ,  then 
the  discrete  maximum  principle  can  be  shown  to  be  equivalent  to 
the  multi-level  system  theory.   However,  the  multi-level  system 
.aory  employs  a  price-adjustment  rule  for  adjusting  the  multi- 
pliers to  achieve  the  optimum  of  the  original  integrated  prob- 
lem.  It  also  gives  an  economic  interpretation  to  the  subprob- 
lems resulting  from  the  decomposition.   But  the  discrete  maximum 
principle  neither  employs  any  method  for  adjusting  the  adjoint 
variables  to  achieve  the  optimum  of  the  problem  nor  makes  use 
of  any  economic  interpretation  to  the  decomposition. 
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In  proving  that  the  adjoint  variables  and  the  Lagrange 
multipliers  are  the  same,  considering  only  the  simple  feedback 
discrete  cases,  there  will  be  little  loss  of  generality.   The 
simple  feedback  discrete  cases  have  been  treated  in  detail  in 
Chapter  III,  Part  One,  by  the  multi-level  system  theory  and  by 
the  maximum  principle  in.  section  1,  Chapter  II,  Part  Two. 

By  comparing  the  necessary  conditions  derived  from  the  dis- 
crete maximum  principle,  equations  (2.2Lj_)  through  (2.29),  with 
those  derived  from  the  multi-level  system  theory,  equations 
(1.10)  through  (1.15)  in  Part  One,  we  see  that  the  Lagrange 
multipliers,  (or  the  prices  of  xn) ,  pn,  in  the  multi-level 
system  theory  are  actually  the  adjoint  variables  zn  in  the  dis- 
crete maximum  principle,  and  that  the  adjoint  function,  equa- 
tion (2.2i|_),  of  the  maximum  principle  Is  actually  a  necessary 
condition  that  the  Lagrangian.  be  stationary  with  respect  to 
its  arguments  xn   ,  equation  (l.ll)  in  Part  One. 

By  comparing  the  computational  schemes,  It  is  easy  to 
recognize  that  the  price-adjustment  rule,  equation  (1.17)  in 
Part  One,  used  in.  the  second-level  problem  of  the  multi-level 
theory  Is  exactly  the  same  as  the  rule  for  adjusting  the  ad- 
joint variables,  equation  (2.21),  suggested  in  the  preceding 
sections . 

In  the  feasible  decomposition  method  of  the  multi-level 
system  theory,  derived  in  section  3,    Chapter  V,  Part  One,  it 
is  assumed  that  the  decision  vector  0n  has  at  least  as  many 
components  as  the  state  vector  xn  and  the  Jacoblan  matrix 
(fn)  n  is  of  full  rank  of  all  6n  and  xn.   It  is  easy  to 
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recognize  that  if  the  components  of  vector  0n  and  x  are  exactly 
the  same,  the  feasible  decomposition  method  and  the  steepest 
ascent  of  the  Hamiltonian  method  suggested  in  the  preceding 
section  will  be  identical.   However,  if  the  components  of  the 
decision  vector  0n  are  fewer  than  the  components  of  the  state 
vector  xn,  then  the  feasible  decomposition  method  fails.   This 
appears  to  be  the  weak  point  of  the  feasible  decomposition 
method.   On  the  other  hand,  if  the  components  of  the  decision 
vector  9n  are  more  numerous  than  the  components  of  the  state 
vector  xn,  the  steepest  ascent  of  the  Hamiltonian  method  fails. 
This  means  that  the  system  will  be  overdetermined  when  all  the 
decision  variables  are  specified,  that  is,  the  number  of  the 
decision  variables  will  be  greater  than  the  number  of  the  in- 
dependent variables.   In  this  case  we  should  use  the  feasible 
decomposition  method.   Or,  if  the  state  variables  are  inter- 
changed with  the  decision  variables,  the  method  of  the  steepest 
ascent  in  the  Hamiltonian  space  can  be  used. 

So  far  we  have  compared  only  the  two-level  structure  of 
the  multi-level  system  theory  with  the  discrete  maximum  princi- 
ple for  the  simple  feedback  process.   It  is  plausible  to  extend 
this  comparison  and  identification  to  the  multi-level  structure, 
which  is  more  complex  than  the  two-level  structure.   However, 
such  an  extension  appears  to  be,  if  not  impossible,  extremely 
difficult.   And  the  use  of  the  multi-level  system  theory  alone 
to  optimize  the  complex  multi-level  structure  also  appears  to 
be  very  tedious  if  compared  with  the  use  of  the  maximum 
principle . 
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Although  the  multi-level  system  theory  and  the  discrete 
maximum  principle  can  be  proved  to  be  mathematically  equivalent, 
the  way  of  approach  to  an  optimization  problem  according  to 
each  method  is  quite  different.   By  means  of  the  multi-level 
system  theory,  we  can  decompose  the  whole  integrated  system 
into  subsystems  with  smaller  dimensions.   The  subsystems  then 
may  be  solved  by  any  of  the  existing  optimization  techniques. 
Thus  it  is  plausible  to  construct,  by  combining  the  multi-level 
approach  with  the  discrete  maximum  principle,  a  powerful  method 
for  optimizing  the  highly  dimensional  complex  system.   This 
should  be  an  area  for  future  work. 


PART  THREE 


PROCESS  ANALYSIS  AND  DESIGN  OP  A  SEQUENTIAL 
REVERSE- OSMOSIS  WATER  PURIFICATION  SYSTEM 
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CHAPTER  I.   INTRODUCTION 

The  reverse-osmosis  water  purification  process  consists  of 
raising  the  pressure  of  an  aqueous  solution  to  a  pressure  above 
its  osmotic  pressure  and  bringing  it  into  contact  with  a  selec- 
tive membrane  which  is  much  more  permeable  to  water  than  to  the 
impurities  (solutes). 

The  use  of  reverse  osmosis  for  water  purification  is  being 
considered  for  saline  water,  brackish  water,  and  process  waste 
xvater.   The  process  is  also  being  considered  for  water  purifi- 
cation in  remote  areas  where  only  small  quantities  of  water  are 
needed.   The  process  analysis  and  design  study  in  this  paper 
are  primarily  intended  for  those  applications  where  saline  and 
brackish  waters  are  to  be  purified. 

Because  of  its  simplicity  and  low  energy  requirements  (the 
water  undergoes  no  phase  change  and  temperature  changes  are 
small),  reverse  osmosis  has  been  drawing  widespread  favorable 
attention  as  a  purification  technique.   It  is  now  well  estab- 
lished that  synthetic  osmotic  membranes,  made  of  cellulose  ace- 
tate, formamide,  and  acetone,  can  be  produced  which  are  highly 
permeable  to  water  and  sufficiently  impermeable  to  dissolved 
salts.   Although  the  reproducibility  and  durability  of  these 
membranes  are  still  in  doubt,  the  results  obtained  to  date  are 
sufficiently  encouraging  to  warrant  a  closer  look  at  the  pos- 
sible economics  of  reverse  osmosis  as  a  water-desalting  process. 

This  study  is  an  attempt  to  investigate  the  reverse-osmosis 
process  in  order  to  find  ways  in  which  the  design  can  be  improved 


'low  model  of  the  reverse-osmosis  unit  is  d  - 
vised.   It  is  based  on  boundary-layer  theory  and  one-dimensional 
.ion  I     ".   A  set  of  system  equations  or  a  system  model 
jh  relates  the  flow  rate,  en       .  quirement,  and  cost  is 

It  is  often  desirable  to  formulate  a  model 
lates  the  behavior  of  the  process.   This  model  may 
be  used  to  find  the  optimum  design  and  operating  conditions 
1  the  system.   The  design  of  a  system  of  reverse  osmosis  units 
connected  in  simple  sequence  is  investigated  in  this  study. 

The  proposed  sequential  process  is  described  and  iLS  advan- 
tages are  c     -Datively  discussed  in  section  2.   In  section  3 
a  boundary-layer  flow  model  is  derived  and  the  process  based  on 
bhis  model  is  analyzed  for  calculating  the  flow  rate  of  fresh 
water  through  the  membrane.   In  section  Lj_,  the  power  requirement 
for  the  process  Is  determined  and  the  cost  function  Is  derived. 
A  conceptual  design  of  the  process  is  given  in  section  $. 

CHAPTER  II.   DESCRIPTION  OP  PROCESS 

The  simplicity  of  the  reverse-osmosis  process  is  apparent 
upon  inspection  of  a  flowsheet  of  the  Aerojet  pilot  plant  (33). 
Basically  the  process  consists  of  a  pumping  system  to  raise  the 
pressure  of  the  brine  solution,  and  of  an  array  of  selective 
membranes.   The  only  energy  consumption  required  by  the  process 
is  that  for  driving  the  pumps.   A  reduction  in  energy  consump- 
tion will  reduce  the  cost  of  the  fresh  water  produced.   The 
minimum  energy  requirement  in  an  ideal  reverse-osmosis  process 
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would  be  achieved  by  applying  a  differential  pressure,  A~P , 
across  the  membrane.   In  other  words,  the  pressure  difference 
should  be  only  infinite simally  greater  than  the  osmotic  pressure 
of  the  brine  solution.   The  concentration  of  the  brine  solution 
should  be  allowed  to  increase  only  inf initesimally  in  the  pro- 
cess.  A  blow down  turbine  can  be  used  to  recover  the  energy  of 
the  high  pressure  reject  brine  solution.   However,  in  the  real 
process  there  is  an  energy  loss  due  to  increasing  the  pressure 
of  the  main  brine  solution,  above  the  osmotic  pressure  and  then 
rejecting  it  to  atmospheric  pressure.   Therefore  there  is  a  min- 
imum energy  requirement  for  the  process,  which  is  different  from 
that  for  the  ideal  process.   Furthermore,  the  fresh  water  flux 
through  the  osmotic  membrane  is  a  function  of  the  pressure  dif- 
ference across  it,  which  is  the  so-called  driving  force.   To 
minimize  the  capital  cost  of  the  separating ~unit  requires 
separating  pressures  substantially  above  the  osmotic  pressure 
of  the  most  concentrated  brine  in  the  system.   Thus  an  economic 
balance  between  energy  and  capital  costs  must  be  achieved  by  the 
proper  selection,  of  operating  pressure  and  reject-brine  concen- 
tration. 

As  we  have  just  stated,  the  pumping  work  is  related  to  the 
osmotic  pressure  of  the  brine  solution.  The  energy  consumption 
is  proportional  to  the  pressure  required  in  each  stage.  The 
higher  the  osmotic  pressure  in  the  separator  units,  the  greater 
the ' energy  consumption.  The  osmotic  pressure  of  sea  water  con- 
taining 35^000  p. p.m.  (i.e.,  3.5  wt.  per  cent)  total  salts  is 
approximately  2l\.   atm.  (Lj-7)  •   In  diluted  sea  water  it  is  roughly 
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proportional  to  the  salt  concentration.   As  fresh  water  is  re- 
moved from  the  brine  solution  the  salt  concentration  at  the 
membrane  boundary  becomes  higher  than  that  of  the  bulk  solu- 
tion.  Accordingly  the  effective  osmotic  pressure,  which  is  the 
osmotic  of  the  solution  at  the  membrane  surface,  increases. 
This  effect  of  salt  build-up  at  the  membrane  surface  is  signif- 
icant with  present-day  membranes.   It  may  become  an  increasingly 
more  important  problem  as  better  membranes  are  developed.   This 
boundary-layer  effect  will  be  discussed  in  more  detail  later. 
To  reduce  this  boundary-layer  effect,  the  concentrated  brine 
solution  in  the  boundary  layer  should  be  mixed  with  the  main 
stream  at  frequent  intervals  or  the  bulk  solution  flow  rate 
should  be  increased  in  order  to  reduce  the  thickness  of  the 
boundary  layer.   This  may  be  accomplished  by  using  a  recycle 
flow  to  keep  the  flow  in  a  turbulent  condition. 

In  an  attempt  to  optimize  the  design  of  a  reverse-osmosis 
process,  Lonsdale,  et  aJL.  {3k-)    considered  a  simple  one-separator 
unit  operation.   The  discharge  salt  solution  in  their  best  re- 
sult contains  nearly'  6  wt  per  cent  of  salt.   To  reduce  the 
boundary-layer  effect, they  proposed  the  use  of  a  recycle  flow. 
However,  from  the  viewpoint  of  thermodynamics,  it  is  unwise  to 
mix  6  wt  per  cent  solution  with  3  •  f?  wt  per  cent  (average)  sea 
water.   Thus  we  propose  a  multi-stage  sequential  system  as 
shown  in  Pig.  1. 

Since  the  osmotic  pressure  is  proportional  to  the  brine 
concentration  in  diluted  brine  solutions,  it  is  advantageous  to 
use  a  low  pressure  in  the  first  stage  where  the  concentration 
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of  salt  is  the  lowest.   Tho  brine  concentration  increases  as 
the  brine  solution  progresses  throu,     e  plant.   Therefore  we 
shall  use  a  stepwise  increase  in  pressure  from  stage  to  s; 
that  is,  the  plant  will  operate  at  as  nearly  an  ideal  reverse- 
osmosis  process  as  possible.   Thus  we  insert  a  high  pressure 
pump  between  stages,  which  is  used  to  increase  the  pressure  in 
ie  process  from  one  stage  to  the  next. 

Instead  of  mixing  the  high  concentration  output  with  the 
relatively  lower  concentration  inlet,  we  propose  to  use  recycle 
flow  at  each  stage.   Thus  we  insert  a  bring  circulation  pump  to 
recycle  the  flow  in  each  stage. 

For  each  stage,  we  use  a  conceptual  shell-and-tube  design 
arrangement  which  clearly  gives  a  lower  capital  cost  per  unit 
membrane  area  than  does  the  plate-and-f rame  configuration. 

CHAPTER  III.   ANALYSIS  OF  PROCESS 

1.   Boundary-layer  Flow  Model 

Water  passing  through  the  membrane  is  supplied  to  the  mem- 
brane boundary  by  bulk  flow  of  solution  normal  to  the  membrane. 
Salt  is  carried  along  with  the  water.   If  a  steady  state  is  to 
be  maintained  without  an  accumulation  of  the  salt  on  the  mem- 
brane, this  salt  must  diffuse  back  into  the  main  bulk  solution. 
A  salt  concentration  gradient  is  established  near  the  membrane 
boundary  such  that  the  net  salt  flux  through  the  membrane  is 
zero.   This  means  that  the  effective  osmotic  pressure  is  greatei 
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than  that  of  the  bulk  solution.   The  situation  near  the  mem- 
brane boundary  is  shown,  in  Fig.  2.   If  we  insert  a  circulation 
pump  at  each  stage,  it  is  possible  to  make  the  flow  in  the  tubes 
be  in  the  region  of  fully  developed  turbulent  flow,  that  is,  the 
bulk  solution  flowing  parallel  to  the  osmotic  membrane  surface 
is  well  mixed,  and  the  existence  of  concentration  and  velocity 
gradients  is  mainly  confined  to  the  laminar  boundary  layer. 

In  the  absence  of  chemical  reaction  the  ratio  between  the 
concentration  boundary-layer  thickness,  6C,  and  the  momentum 
transport  boundary-layer  thickness,  6,  is  shown  by  Bird,  _et  al. 
(35)  to  be  a  constant  which  is  dependent  only  on  the  value  of 
the  Schmidt  number,  i.e., 

6 
or 


C  =  Sc"1/3 


6C  =  8  •  Sc--7J  (1) 


rV3 

where  Sc  is  the  Schmidt  number  which  is  equal  to  the  kinematic 

viscosity   ,  divided  by  the  diffusion  coefficient  for  the 
solution,  DQ,  or 

Sc  =  —  .  (2) 

The  thickness  of  the  laminar  sublayer  for  turbulent  flow 
through  pipes  is  given  in  Schlichting  (l\.)    to  be  equal  to 

5-5—  (3) 

v-::- 

sq  cm 

where  ))   is  the  kinematic  viscosity  ( )  ,  and  v#   is  the 

sec 

friction  velocity  defined  as 
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Fig.  2.   Velocity    and    salt  concentration' gradients 
in    boundary    layer   adjacent  to  a  membrane. 
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in  which  t^  is  the  shearing  stress  at  the  wall.   v_;;_  is  a  measure 
of  the  intensity  of  turbulent  eddying  and  of  the  transfer  of 
momentum  due  to  these  fluctuations.   For  turbulent  flow  through 
pipes,  the  shearing  stress,  Tq ,  can  be  calculated  from  the  fol- 
lowing equation  (I4.)  : 

t0  =  0 .0225  /V 'A  (i)i/k  (5) 

R 


where  f   is  the  density  of  the  fluid  and  U  the  maximum  velocity 

which  is  proportional  to  the  mean  velocity  u.   The  ratio  ~ 

U 
equals  0.8  for  turbulent  flow  through  pipes.   R  is  the  radius 

of  a  pipe  and  d  =  2R  denotes  the  diameter  of  the  pipe.   Intro- 
ducing d  into  equation  (5>)  and  then  substituting  equations  (I4.) 
and  (5)  into  equation  (3)  yields 
5  D  $)) 


s  = 


S  /0>0225u7/^   (JL)1A 

"  d    2 


5 


-tf^IiF  .    2i/8  .    (  —  )7/8  .   d-l/8  .   v<i/8-D 

0.8 

5   •    (0.8)?/8  d 


0.15   •    L09P  d?/8  )T7/ti 
33.3    •    (0.82^)    d 

1.09    (^j-)7/8 


25.2   d  25.2   d 


^d~^_Rj7^ 
( ) 


(6) 


M- 
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ud   ud 
where  Re  =  —  =  is  the  Reynolds  number.   Substituting 

equation  (6)  into  equation  (l),  we  obtain 
25.2  d 
°C  =  Scl/3  Re7/8  '  (?) 

In  a  sequence  of  stages  the  concentration  and  velocity 
differ  from  stage  to  stage .   If  the  viscosity,  density,  and  dif- 
fusivity  are  assumed  to  be  constant,  and  the  diameter  of  the 
tubes  is  assumed  to  be  the  same  in  each  stage,  then  the  velocity 
will  change  at  each  stage.   Thus  the  Reynolds  number  will  be 
different  in  each  stage.   The  concentration  boundary-layer 
thickness  at  the  nth  stage  then  becomes 

25.2  d  -, 

6"  =  —-7- -77;  (7a)1 

C   Sc1^  (Ren)7/8 

From  equation  (7)  it  can  be  recognized  that  the  larger 
Reynolds  number  which  can  be  achieved  by  increasing  the  circu- 
lation rate  reduces  the  boundary-layer  thickness.   Thus  the  in- 
crease in  osmotic  pressure  arising  because  of  the  boundary  layer 
will  undoubtedly  have  to  be  controlled  by  providing  adequate 
circulation  rates  through  the  tubes.   We  have  assumed  that  the 
system  is  isothermal  and  that  there  is  no  precipitation  of  salt 
on  the  membrane  surface.   The  experimental  result  of  Merten 
(37)  suggests  that  alternative  procedures  may  be  devised  to 
control  salt  precipitation.   The  effects  of  circulation  on  flow 
through  an  osmotic  carrier  have  been  experimentally  observed  by 


Note  that  the  superscript  number  n  indicates  the  stage 
number.   Exponents,  where  required,  will  be  written  outside  of 
parentheses  or  brackets. 
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Merten  (37) •   His  results  are  in  good  agreement  with  equation 
(7). 

2.   Simple  Sequential  Multi-stage  System 

The  proposed  sequential  reverse-osmosis  water  desalination 
system  shown  in  Pig.  1  is  of  the  form  of  a  simple  sequential 
multi-stage  model  (39)  as  shown  in  Fig.  3.   Each  stage,  except 
the  last  stage,  includes  a  membrane  separator  unit,  a  recycle 
pump,  and  a  high  pressure  pump  between  stages.   Figure  I4.  gives 
a  schematic  representation  of  the  nth  stage.   The  last  stage, 
stage  N,  includes  a  blowdown  turbine  in  its  outlet. 

Let  us  now  define  the  following  symbols: 

x11  =    the  average  mass  fraction  of  salt  at  the  nth 
stage 

qn  =  the  mass  flow  rate  of  the  brine  solution  dis- 

lb 
m 

charged  from  the  nth  stage  ( ) 

hr 

Wn  =  the  mass  flow  rate  of  fresh  water  product  from 

lb 
m 

the  nth  stage  ( ) 

hr 
W~  =  the  total  mass  flow  rate  of  fresh  water  produced 

lbm 

from  the  whole  system,  ( ),  that  is, 

hr 
N 

wf  =  "ZT  wn  '(8) 

n=l 

N   =  the  total  number  of  stages  in  the  sequence  of 
the  process. 
Here  we  assume  that  the  salt  concentration  of  fresh  water  pro- 
duced is  negligible.   We  also  assume  that  the  recirculation  rate 
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Fig.  4.     The   representation   of  a    single    stage. 

Where     x    represents   the   state    vector . 

0  represents  the  decision   vector. 
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is  sufficiently  high  with  W  .   Therefore  we  may  also  assume 
that  the  salt  concentration  in  the  bulk  solution  is  constant, 
that  is,  xn  =  xn,  where  xn  is  the  salt  concentration  of  the 
brine  stream  leaving  stage  n. 

Then  the  total  material  balance  for  the  process  as  a 
whole  is 

q°  =  qN  +  Wf  .  (9) 

The  salt  material  balance  for  each  stage  is 

„0  0  _  „lvl  _     _  ^nvn  _  „uvn  (m) 

q  x  —  q  x  —  .  .  .  —  qx—  qx.  ^  iu  ; 


3.   The  Volumetric  Plow  Rate  of  Fresh  Water 
Produced  at  the  nth  Stage,  Fn 


The  volumetric  flux  of  water,  Fn,  through  a  membrane  of 
constant  permeability  has  been  reported  by  Merten  (37)  as 

ft3 

Fn  =  K(APn  -  TtS)  ( )  (11) 

ft2  -  hr 

ft-3 

where  K  =  the  membrane  constant  ( ) 

ft2  -  hr  -  psi 

Apn  =  the  pressure  difference  across  the  membrane  of  the 

nth  stage  (psi) 

7t«  =  the  osmotic  pressure  of  the  brine  solution  at  the 

membrane  surface  (psi). 

To  relate  the  osmotic  pressure  to  the  brine  concentration, 

Merten,  _et  al.  (!|.0),  has  found  that  the  expression 

Tig  =  12,100  x^  (psi)  (12) 

fits  the  experimental  data  of  Tribus,  et  al.  (lj.1),  where  x   is 

the  average  mass  fraction  of  salt  concentration  at  the  membrane 
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surface  of  the  nth  stage. 

Now  let  us  apply  the  boundary-layer  flow  analysis  to  the 

nth  stage.   The  pressure  is  uniform  in  the  high  pressure  chamber 

if  the  pressure  drop  due  to  bulk  flow  is  negligible.   A  salt 

material  balance  inside  the  concentration  boundary  layer  for  a 

plane  parallel  to  the  membrane  is  described  as  follows. 

dx11        xn 

-D  =  Fn  (13) 

8  dy       1  -  xn 

jA  n 
dx 

where  Do  is  the  rate  of  migration  of  salt  component  in  the 

8  dy 

direction  from  the  membrane  surface  to  the  main  bulk  solution 

xn 

by  diffusion,  and  Pn  is  the  volumetric  flow  rate  of  salt 

1  -  5cn 

in  the  direction  from  the  main  bulk  solution  to  the  membrane 
surface  by  bulk  motion  or  convection. 
Therefore 

dxn   Fn    xn 


dy    D   1  -  xn 


(34) 


Since 


^n 
xu 

xn<<;l,   =  xn 

1  -  xn 

equation  (ll\.)    can  be  simplified  to 

dxn   Fn 

-  —  =  —  dy  .  (15) 

.  xn    Da 

If  we  assume  D   to  be  independent  of  x  ,  integrating  the  right- 


a 

n 
hand  side  of  equation  (if?)  from  0  to  5 

„    ^n    ,   An   .  _ ■ 
from  x„  to  x  yields 


and  the  left-hand  side 
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An 

r  xs  dxn 


^n 
x 


\ 


n 


pn 


n 


dy 


'0 


or 


xn   Fn 

/n  —  =  —  6   . 

X*    D    C 
a 


(16) 


Substituting  the  concentration  boundary-layer  thickness,  equa> 
tion  (7a),  into  equation  (16),  we  have 

x^   Pn        25.2  d 

x^~~  D~   (Sc)1/3  (Ren)7/6  ' 


£n   = 


(17) 


If  we  use  the  approximation 

An     An 

/n  —  =  ( 1) 

xri    .  £n 


(18) 


n 


in  equation  (17)  and  solve  for  F  ,  we  obtain 


xn       (Sc)1/3  (Ren)7/8  D. 


Fn  =  (_£  .  1} 
x' 


(19) 


25.2  D 

Substituting  equations  (19)  and  (12)  into  equation  (11)  and 

solving  for  Xo,  we  obtain 

25.2  d 

xn  (1  +  KAPn  — -r- ) 

(Sc)1/3  (Ren)7/6  Da 
.       (20) 


Xn   = 


1  +  12,100  xn  K 


25.2  d 


(Sc)1/3  (Ren)7/8  D, 


Substituting  equation  (20)  into  equation  (19),  we  obtain  the 
volumetric  flow  rate  of  fresh  water  per  unit  area  of  membrane 
from  the  nth  stage  as 
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/^n.' 


K(APn  -  12,100  x 

Fn   =  - 

25.2  d 
1  +  12,100  xn  K 


(Sc)1/3  (Ren)?/8  D, 


K(APn  -  12,100  xn)     ft3 

(— )  (21) 


xn       ft2  -  hr 


1  +  C 


(Re*)?/8 

where 

25.2  d 
C  .=  12,100  K 


(Sc)i/3 


D 


a 


K  d 

=  3.0£  x  1(P  ™ .  (22) 

(Sc)1/3  Da 

The  overall  material  balance  around  the  nth  stage  is 

qn  =  qn-l  _  wn  #  (23) 

Thus  the  fresh  water  produced  from  the  nth  stage  is 

Wn  =  Fn  f>   Sn  (23a) 

where  Sn  is  the  membrane  area  of  the  nth  stage  (ft  )  and  f   is 

lbm 
the  density  of  the  fresh  water  ( ) . 

ft3 

The  salt  balance  around  the  nth  stage  is 

0  0    n^n 
q  x  =  q  x 


or 


xXi  = 


q11 


0  0 
q  x 


qn-l  _  Fn  f>   sn 
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q°x° 


x°q° 

Fn  P   Sn 

fcn-1 


X0xn-1 


Sn 
x°  -  x11"1  FnP(—) 
q° 


(21+) 


I4..   Relation  Between  the  Reynolds  Number  Ren 
and  the  Recycle  Ratio  Rn 


The  cross- sectional  area  through  which  the  brine  solution 

passes  at  the  nth  stage,  An  (note  that  this  is  different  from 

the  membrane  area  Sn) ,  depends  on  the  design  of  the  separator 

unit  at  each  stage.   Often  it  is  economical  to  use  unifying 

equipment  at  each  stage;  thus  we  use  at  every  stage  a  set  of 

m-tubes  arranged  parallel  to  each  other  as  a  separator  unit. 

Then  the  cross-sectional  area  remains  the  same  at  each  stage  and 

can  be  given  in  terms  of  the  diameter  of  the  tubes  as 

rrm(d)2 
An   =  A  =  ,  (25) 

k 
Similarly,  Sn  -  S  for  all  stages. 

The  fluid  velocity  inside  the  tubes  of  the  nth  stage  is 
qn-l  (1  +  r") 


An 

u   = 


AP 

^qn_1  (1  +  Rn) 


nm(d)2  f 


(26) 


The  relation  between  the  Reynolds  number  Ren  and  the 


recycle  ratio  Rn  at  the  nth  stage  is 
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du7        d/>    W1"1    (1  +  Rn) 
Ren  _  = 


[L  u-  PmJi(d) 

l^q""1    (1    +   Rn) 

=  .  (27) 

u-rrmd 

n-l         1°X° 

Substituting  q  =  into   equation    (27),    we    obtain 

xn-l 

i  0  0 

Ren  =   . (1  +  Rn)  .  (28) 

\im   xn_1  dit 


CHAPTER  IV.   COST  ANALYSIS 


1.   Operating  Cost 

Energy  requirements  per  lb   of  fresh  water  product  can  be 
determined  as  follows. 

Let  us  consider  any  two  successive  stages  as  shown  in 

n 
Pig.  '5-   E   represents  the  pumping  work  of  the  high  pressure 

n 
pump  at  the  nth  stage  and  Ep  represents  the  pumping  work  of  the 

circulation  pump  at  the  nth  stage. 

Let  Y[    ,  h  ,    and  t\      be  the  mechanical,  pump,  and  turbine 
efficiencies,  and  Y\  -  the  loss  factor. 

A .   Energy  Requirement  for  the  High  Pressure  Pump  at  the 
nth  Stage .   The  pumping  work  E-,  is  primarily  used  to  increase 
the  pressure  from  P     to  Pn.   Since  the  velocity  difference 
between  the  two  successive  stages  is  small,  the  kinetic  energy 
losses  and  friction  energy  losses  can  be  included  in  the  pump 
efficiency.   Thus  the  power  requirement  for  high  pressure 
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pumping  at  the  nth  stage,  E-,  ,  is 

1  +  Y[~      Pn  -  Pn_1     _   psi  -  ft3 

E?  = qn_1  ( )  .  (29) 

1     nm>lV  e  hr 

Fop  one  lb   of  total  fpesh  water  produced,  we  have 
m 

E?   1  +  nr     P"  -  Pn_1   qn_1   P^  -  ft3 

—  = - — t— ( )  .      (30) 

Wf    Wp       r       Wf      lbm 

Substituting  equations  (9)  and  (10)  into  equation  (30)  yields 

0  0 
q  x 


■l 

i  +  Vf 

Pn   -    P*"1             x^1 

wf 

^rrtfp 

f                              x° 

q0(1"^       ■ 

1   +7?f 

nmnP 

pn   _   pn-l                 x0 

e         1      xo 

xn-l(l   _   _) 

XN 

i  +  nf 

>1m>?p 

A?      -  AP*"1                  x° 

f3                                        x0 

xn-l(l ) 

XN 

Enepgy   Requirement    fop    the    Recycle    Pump 

(3D 


Stage .   The  enepgy  required,  Ep,  includes  the  enepgy  of  circu- 
lating the  qn_1Rn  lkm/hr  of  fluid  and  that  of  the  q11"1  flow  wopk. 
The  fpiction  loss  comes  lapgely  from  the  fluid  flowing  in  the 
membrane  separator  unit.   This  lost  wopk  is  (lj.2) 

„  (un)2    l        _  i  +  nf 

E"  =  kl   (-)  qn"1(l  +  Rn) £  (32) 

d  2gc    d  V>?p 

where  f  is  the  friction  factor,  L  is  the  length  of  the  membrane 

separator  unit,  and  d  is  the  tube  diameter  of  the  membrane 

separator  unit.   qn_1(l  +  Rn)  is  the  amount  of  fluid  flowing 
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through  the  membrane  separator  unit,  which  is  equal  to 

qn-l(l  +  Rn}  =  m  1_   'un  (33) 

k 
where  m  is  the  total  number  of  tubes  inside  each  membrane 

separator  unit  as  mentioned  previously. 

As  discussed  before,  the  fluid  flow  within  the  membrane 
separator  chamber  would  be  in  turbulent  flow.   Under  this  con- 
dition, the  friction  factor  can  be  approximated  by  (I4.3) 

0.01+6 

f  =  o-p   •  (3^ 

(Ren)U*^ 

Thus  for  one  lb   of  fresh  water  produced,  the  energy  re- 
quired for  the  recycle  pump  is 

E^   11(0.01+6)    (un)2   L    rrm(d)2  unP      1  +  >?f 


2  _ 


(-) 


Wf    (Ren)0'2     2gc    d        l4.Wf      ^m^p 

0.023  dun/°     .         \jl      _    Lrrmd  P      1   +  r\r 

( )3  ( )3  — ~      .         (35) 


(Ren)0-2     [i  dP  gcWf    >Zm^p 

Note  that  Lmrcd  is  equal  to  the  membrane  surface  area  S,  and 
du"/9 


=  Ren  .   By  using  equations  (9)  and  (10),  equation  (35) 

becomes 

-2-=   0.023  (Ren)2*0  ( )J  .        (36) 

Wf  d^     n,     x°    Scttirrtp 

C.   Energy  Recovery  at  Reject-brine  Turbine .   The  equation 
for  the  energy  recovery  from  depressurizing  the  high  pressure 
brine  solution  will  be  of  the  same  form  as  equation  (30)  which 
gives  the  energy  requirement  to  pressurize  the  brine  solution. 
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Thus  we  have 


E,  PN  -  P°    qN 


=  V*m  (1  -  Vf) 


*p'im  x-    ii'     /) 


x' 


0 


Apj     x^ 


Vim  (1  "  *?f)  — 


x° 


(37) 


APN    x° 

=  %nm   (1-  V   —  -i— -  (37) 

where  P   is  the  discharge  pressure  (usually  it  is  one  atmos- 
phere) . 

If  the  energy  required  is  supplied  from  electricity  the 
electrical  power  cost,  Ce,  is  assumed  to  be  $0,005  per  kw-hr 
in  all  cases. 


2.   Capital  Cost 

A.  Pump  and  Turbine  Installation  Cost,  Cp.   For  simplicity 
the  costs  of  pump,  turbine,  and  motor  are  assumed  to  be  directly 
proportional  to  horsepower  rating  in  the  horsepower  range  of 
interest.   An  f.o.b.  cost  of  $100  per  kw  has  been  assumed  by 
Merten,  et  al.  (3ij.)  . 

B.  Membrane  Separator  Unit  Cost .   Because  of  a  lack  of 

information  about  the  cost  of  this  type  of  equipment,  the  cost 

equation  which  Merten, '  e_t  al.  (3^4-)  derived  is  used.   It  is 

W*:         Pmd       AVn                       0.^       0.189      f^~ 
°-=   (— )    (— )    (1.62  +  —  +  —  f  -J3L- )  (38) 


Wn         crm  PF*  l/d  L/D      '    APn 
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where 


Wo  =  the  mass  of  the  shell-and-tube  membrane  separator 

unit  of  nth  stage  (lb) 

lbm 
fm  =    the  density  of  the  material  of  construction  ( -) 


ft3 


cr  =  the  allowable  stress  of  the  material  of  con- 
vm 


struction  (psi) 
L /d  =  the  overall  length-to-diameter  ratio  of  the 
membrane  separator  unit 
Wn  =  the  fresh  water  produced  from  the  nth  stage. 
Changing  equation  (38)  into  the  cost  per  lbm  of  the  fresh  water 
produced,  we  have 


wS      wn      w^        ^md     APn  o.$k      0.189    /en,    wn 

-S  =  —    .   -§  =    ( )( )(1.62    +  — —   +  —  Y  — )    —       . 

Wf        Wf        Wn  crm        Fnf  L/D  L/D  pn      Wf 

(39) 

By  using   equations    (9),     (10),    and    (23a),    equation    (39)    becomes 


Wo         /*md      APn                   0.5J+       0.189     /  crm            FnS^ 
—  -    ( )( Ml.62   +  — —  +  —  J )    - 

Wf  °~m        FUf  L/D  L/D  APn 


XN 


Pmd        s^pn  0.5U-      0.189    /  o-m 

=  - —  (1.62  +  _.  +  _ _ y-2-)   .    (ko) 

crm  jP  l/d  l/d    1    4p* 

q  "    XN 

We  assume  that  the  cost  of  the  membrane  separator  unit,  which 
includes  fabrication  and  installation  costs,  is  proportional  to 
the  weight  of  the  material  used.   And  the  unit  cost  of  the 
material  of  construction  is  Cs  $/lb. 

The  annual  capitalization  charge  for  these  equipment  items 
is  calculated  at  0.0 71+  of  the  initial  cost  per  year,  as 
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recommended  in  the  Office  of  Saline  Water  Report  (ijij.)  .   An 
assumption  of  a  load  factor  of  330-on-stream  days  per  year  gives 
a  capitalization  charge,  *-p  ,  of  9.1\.   x  10    of  initial  cost  per 
hour  on  stream. 

The  total  cost  contributions  of  the  system  are  in  the  form 

n  n 

N        E-i  N        E9 

ct=  (rcp  +  ce)(E    —  +  Z    -> 


n=l     W^        n=l     Wf 


E- 


rn 


J3  #-      W^ 

+    (y.C  Ce)    —  +    ^Cs    2-       -2      .  (1H) 

F  Wf  n=l     Wf 

Substituting   equations    (31),     (3&),    (37),    and    ( I4.O )    into    equation 
(1^.])  and   combining  all   constants,    we   have 
x°  N      APn   -  Apn-1 


Ct   =   Bl  n     ?" 


x°      1  xn-! 

1 

XN 

S     •        xN  .  N  P    o  N. 

+    (— )(— -)    (b2     I      (Ren)2'°   +   B.     I    APn 

q0      XN   _    x0      I    2      !  J      x 

1   +   V?f 
where      Bx  =    (  f  Cp   +  Ce)    (ll3) 

^>V7p 


1  +  >if     ^  ^    p.    ^  3 
Mm' 


b2  =  0.023  (fcp  +  ce)  P{—t;)  (^} 

Mmfp  d  P 


B.    =  —    (1.62   +  — - -)  (l&) 

a~~  L/D 


m 


0.189^CS  Pm  d 
Bh    =  —7= ; (1|6) 

B5   =    (f  Cp   -    Ce)    Ylp  Y\m    (1    -   Ylf)  jf>  .  (kl) 
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CHAPTER  V.   CONCEPTUAL  DESIGN 
1.   Assumptions 

(a)  The  flow  model  is  fixed,  but  the  geometric  parameters 
of  the  membrane  separator  unit  itself,  i.e.,  the  pipe 
diameter  d,  length-to-diameter  ratio  L  D,  and  the 
total  tubes  used  m,  are  chosen  arbitrarily. 

(b)  Each  stage  is  geometrically  similar.   The  membrane 
of  each  stage,  S,  is  identical. 

(c)  Salt  concentration  of  fresh  water  produced  is  assumed 

to  be  equal  to  zero,  i.e.,  salt  components  cannot 
pass  through  the  membrane. 

(d)  Peed  saline  water  solution  contains  3 • 5  weight  per 
cent  of  salt  components  (average). 

(e)  The  system  is  isothermal,  and  there  is  no  precipi- 
tation of  salts  on  the  membrane  surface. 

(f)  The  costs  of  the  pump,  turbine,  and  motor  are  assumed 
to  be  directly  proportional  to  horsepower  rating  in 
the  horsepower  range  of  interest. 

(g)  The  cost  of  the  membrane  separator  unit  is  propor- 
tional to  the  weight  of  the  material  used. 

(h)   D  ,  n,  P   are  assumed  to  be  constant  in  the  concen- 
tration range  of  interest. 

(i)   Membrane  constant  K  is  assumed  to  be  independent  of 
pressure. 


2.   System  Variables  and  Number  of  System 
Variables  (for  the  Total  N  Stages; 
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(a)   Mass  flow  rate  of  brine  solution. 


0    12 
Q  >    Q    >    0.*  •  •  •  t    m. 


qN  ;     N  +  1 


(b)  Salt  concentration. 

-A.   j,   A.   j(    •  «  •  ^   A     j  1M 

(c)  Operating  pressure  at  each  stage. 

N 

(d)  Fresh  water  produced  from  each  stage. 

N 


P1,  p2,  .  .  .,  pN  ; 


W1,  W2,  ..  .,  ¥N  ; 

(e)  Recycle  rate  of  each  stage. 

R1,  R2,  ...,  RN  ;  N 

(f)  Membrane  area  of  one  unit. 

S  1 

Total  number  of  system  variables  =  5>N  +  2 


3.   Relations  and  Number  of  Relations  Among 
Various  Variables 


(a) 


Total  material  balance  at  each  stage. 
.0  _  T.rl  ^   l 


q-  =  w^  +  q 
q1  =  W2  +  q 


N 


qN-l  =  wN  +  qN/   • 
(b)   Salt  material  balance  at  each  stage 


0  0     11 
q  x  =  q  x 


2  2 

q  x   = 


N  N 
=  q  x   ; 


N 
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(c)   Mass  transfer  at  each  stage. 

Sf>  KlAp1  -  12,100  xi; 

W1  =  s^F1  =  


x1 


1  +  C 


W2  =  sfF2    = 


(Rel)^/8 
Sf  K(A  P2  -  12,100  x2) 


x' 


1  +  C 


(ReO 


2)7/8 


J 


N< 


WN  =  S^N  = 


SpK(A  P   -  12,100  xiN) 


x 


N 


1  +  C 


(ReN) 


TJa 


N 


The  total  number  of  relations  is  thus  3N.   The  number  of  inde- 
pendent variables  is  calculated  as  follows: 
Total  variables  -  total  relations 
=  5N  +  2  -  3N  =  2N  +  2  . 
Note  that  from  the  design  equations  and  cost  function  we 
can  see  that  S  and  q°  always  appear  together  in  the  term 
(S/q°).   If  we  use  {S  Icp)    as  a  single  variable  in  place  of  S  and 
q  ,  we  can  reduce  by  one  variable  the  total  number  of  independent 
variables.   Thus  the  total  of  independent  variables  becomes 
2N  +  1.   The  mechanical,  pump,  and  turbine  efficiencies,  Y\    , 
yi-,   Y[Vf    are  all  treated  as  parameters. 
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14-.   Illustrative  Examples 

The  following  data  are  used: 

Da   =  1.5  x  10  ^  cm^/sec,  the  diffusion  constant  for 
Nacl  in  water  at  room  temperature  (I4.6). 

J       -ft3 
K   =  0.86  x  10 -^  (324.) 

ft^  -  hr  -  psi 

V>   =  8  x  10"3  for  the  brines  at  30°C   (I4.7) 

sec 


/>   =  62.24.  lbm/ft3 


1 
d     —  ft 
214- 

L/D  =  8  (3k) 

P  =  2.7  x  62.24.  lbm  ft^  density  of  aluminum   (JL|3) 

(7L  =  15,000  psi,  allowable  stress  of  aluminum   (I4.3) 

(f  =  9.2+.  x  10-6  hr"1 

Cp  =  100  $  per  kw   (32;) 

Ce  =  0.005  $  per  kw-hr   (324-) 

Cg  =  L\.tl\.   $  per  lb   of  aluminum   (324-) 

>?m  =  0-9 

nv   =  \  =  °  • 8 

Y\f      =   0.1  . 
Calculation  of  the  constants  is  carried  out  as  follows: 

-  Scl/3  =  ,JL,l/3  .  (fl3C"'3W3  .  8.x 

Da         1.5  x  10-5 
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C   = 


3.05  x  10-5  Kd 
ScV3    D 


3.05   x   10  £   psi    x  0.86    x   10"^ 


ft- 


ft 


ft2-hr-psi        21+ 


8.1   x   1.5  x  10-£ 


cm2        36OO    sec 


sec 


x 


x  0.001076 


hr 


ft' 


cm 


=    2.32    x   10 


3 


1    +  Vj 


f 


Bl  =    CCp  +  Ce) 

1  $ 

=    (9.1+  x  10"6  —  x  100  — 
hr  lew 


+  0.005 


$ 


kw-hr 


•) 


1   +  0.1 


lb 


in 


x  3.766    x   10-7 


62.1+  x  0.9   x  0.8 

f  t3 

kw-hr  12    in 


=  O.7887   x  10-8    ( 


ft-lb 
$ 


x    (• 


) 


ft 


psi-lb 


in 


-6 


B2   =   0.023    (9.1|   x  10"°    x   100    +   0.005) 


$ 


kw-hr 


1   +  0.1 


x 


lb.p-sec 
(— ) 


0.9   x  0.8  x  32.2      lb™-    ft 


8  x  10-3 


lb 


rn 


cm' 


sec 


x   62. [j.   x    ( 


ft3  1/2I+  ft 


ft* 

x  O.OOIO76  ) 

,2 


3 


cm' 


kw-hr                      sec 
x  3.766   x  10"'    x  3600  


lb^-ft 


hr 
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$ 
=  0.14-835  x  10"17  ( ) 

ft2-hr 

Bo  =   b  m        (1.62  +  ) 

J  <r  l/d 

A     -i         $               lbm    1 
9.1+  x  10"°  hr"1  x  I4..I4.  x  2.7  x  62.1;  x  —  ft 

lbm  ft3   21+ 

15,000  psi 

0.51+ 

(1.62  +  ) 

8 

$ 

=  0.3266  x  10"7  ( ) 

psi-f t  -hr 

0.189  ifCs  rm   d 

Y(Tm     L/D 

,         $  lb    1 

0.189  x  9.1+  x  10-b  x  L|_. JL|_  x  2.7  x  62.1|.  —  x  —  ft 

hr-lbm  ft3    2I4. 

f/l5,000  psi  x  8 

=  0.133  x  10"6  (—= ) 

y  psi-ft2-hr 

B5  =  (rcp  -  ce)  y\vnm  ci  -  W^ 

6  $ 

=  ( 9.1+-  x  10"°  x  100  -  0.005)  x  0.9  x  0.8  x  (1  -  0.1) 

kw-hr 

1  ft3  kw-hr     12  in 

x  x  3.766  x  10"'  x  ( )d 

62.14.  l^m  lbf-ft     ft 

=  -0.2286  x  10-8  ( ). 

lbm  -  psi 

In  this  illustrative  example,  the  average  salt  concentration 

within  the  membrane  separator  chamber  of  the  nth  stage  xn  is 

assumed  to  be  equal  to  the  salt  concentration  in  the  outlet  brine 


914- 


solution  of  the  nth  stages  xn. 

From  equation  (25),  we  see  that  this  assumption  is  valid 

for  high  recycle  ratio;  thus  equation  (21)  becomes 

K(Apn  -  12,100  xn) 
Fn  =  .  .         (50) 

1  +  C  77ft 

(Ren)?/8 

Substituting  equation  (50)  into  equation  (2I|_)  and  simplifying, 

we  obtain 

1     s       K(Apn  -  12,100  xn) 


-n 


xn-l   q0   x0  xn 

1  +  C  - 


=   i  .     (5D 


(Re^)?/8 

A.   Calculation  of  the  Cost  for  ja  3-stage  System  for 
Various  Recycle  Ratios  but  with  the  Same  Operating  Pressure  at 
Each  Stage . 

The  total  number  of  independent  variables  is 

2N  +  1  =  7  • 
By  choosing  the  Reynolds  number  for  each  stage  as 
Re1  =  .97  x  10^ 
Re2  =  1.9Z+.  x  10^ 
Re3  =  2.9  x  10^  , 
the  operating  pressure  drop  as 

P1  =  p2  =  p3  =  loll;.  7  p3i 
or     AP1  =  4P2  =  Ap3  =  1,000  psi, 
and  the  discharged  concentration  as 

x3  =  0.071, 
then  the  whole  system  will  be  fixed.   Note  that  using  the  Rey- 
nolds number  at  each  stage  as  an  independent  variable  renders 
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the  calculations  easier  than  using  the  recycle  ratio  as  an  in- 
dependent variable. 

The  values  of  the  Reynolds  number  used  are  based  on  the 
assumption  that  the  fluid  velocities  in  the  separator  chamber 
are 

U1  =  2  ft/sec 

U2  =  ij.  ft/sec 

Ij3  =  6  ft/sec. 
The  corresponding  values  of  (Ren)  ''  ,    n  =  1,  2,  3 ,    are 

(Re1)7/8  =  (0.97  x  10^) 7/8  -  3.08  x  103 

(Re2)7/8  =  (1.91).  x  10^) ?/8  =  5.61  x  103 

(Re3)7/8  =  (2.9  x  101*-)7/8  =  8.0X4.  x  103  . 
Substituting  the  known  values  into  equation  (5l),  we  have 

lbY 

0    .  62. k 

S  ft2 

(  — ) 


'm 


x 


ft3 


0.035  q°      lt>m/hr 

6     ft3 
0.86  x  10_b 


f t2-hr-psi 


0.035 


(1000  -  12,100  x1)psi 


=   1 


1   +   2.32   x  10 


3 


X 


1 


3.08  x  10 3 


and 


x' 


ft' 


62.  k 


lb 


m 


ft 


1  S 

txl  q°      lbm/hr        0.035 

.  ft3 

0.86  x  10  "^ 


3 


f t2-hr-psi 


1  +   2.32   x  10 


(1000    -    12,100    x2)psi 


=   1 


3 


xc 


5.6l  x  103 
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and 


(1  S     ft 

0.071 ( — ) 


62.14. 


lb. 


ft 


x2     q°   lbm/hr    0.035 

ft3 

(1000  -  12,000  x  0.07D 

ft2-hr-psi 


0.86  x  10"^ 


=  1  . 


1  +  2.32  x  10 


-3 


0.071 


9.0[[  x  103 

By  trial  and  error  we  can  solve  these  three  equations  with  three 

unknowns.   The  results  are 

S  ft2 

—  =  0.11 

nO 


lbm/hr 


.1  _ 


x-  =  0  .OI4.7 

x2  =  0.06  . 
Then  we  can  calculate  the  cost  for  this  system.   First  we 
calculate 


Z.     (Ren)2'8  =  (0.97  x  10^)2,8  +  (1.914  x  106)2'8 
+  (2.9  x  10^)2'8  =  I4..385  x  1012 


n=l 


£"  (ZiPn)  =  3,000  psi 
n=l 


£-  (AP")1/2  =  3  x  (lOOO)1'2  =  91+..92  l/p~sl  . 


n=l 


Substituting  these  values  into  equation  ( l_j_2 )  ,  we  have 


Ct  =  Bl 


X°    N   APn  -APn_1     S 

t  : +  <-r) 


,N 


xO  1 

1  -  — 

XN 


rn-l 


q°   xN  -  x° 


/    N 

(b2  T    ( 


Re 


m  2.8 
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0 


N  N  -,  /o)  xL 

+  Bo    T.    APn  +  B,      f     (APn)V2 /+  B.  ApN 

3      1  ^      1  j  ^   XN    _    x0 

o  $  0.035  1000    psi 

O.7887   x  10"°  ( +   0    +   0) 

psi    -    lbm  0.035        0.035 

1  -  

0.071 

ft2        0.071  ir,        $ 

+  0.11  ( )(0.1j.83 5  x  10-1' 


lb^/hr  0.071  -  0.035)  ft2-hr 


$ 

x  I4.. 385  x  1012  +  0.3266  x  10"7  ■  x  3000  psi 

psi-f t2-hr 

*  /— 

+  0.133  x  10"b  -== x  914-.92  /psi 

ypsi-ft2-hr 

«    $  0.035 

-  0.2286  x  10"°  x  1000  psi  x 


lb  -psi  0.071  x  0.035 


m 


=  14-.2255  x  10"£  =  0.35^5  ( )   . 

lbm  1000  gal 

B.   Calculation  of  the  Cost  for  a_  System  of  Three  Stages 
with  a_  Stepwise  Increase  in  the  Recycle  Ratio  and  Operating 
Pressure  from  Stage  to  Stage ♦ 

As  in  (a),  the  system  has  seven  independent  variables. 
Here  we  use  the  same  values  of  the  Reynolds  number  as  in.  (a). 
However,  the  operating  pressure  is  increased  from  stage  to  stage, 
as  follows. 

AP1  =  1000  psi 

AP2  =  1250  psi 

Ap3  =  1500  psi. 
The  discharging  concentration  is 

x3  =  0.07  . 
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Then    the    system    is    fixed.       Substituting   these    values    into    equa- 
tion   ( 5l) ,    we   have 


f      1 


S        62. If        0.86    x   10"^    (1000    -    12,100    x1) 


-    (  — ) 


=   1 


0.035  q°      0.035 


xJ 


1   +    2.32    x   10- 


3.08   x   103 


1             S        62. If        0.86    x   10"^    (1250    -    12,100    x2) 
x2\  —   -    (  — )    f 


i1  q°      0.035 


1   +    2.32    x   10 3 


5.6l   x   10 3 


0.07 


fl  S        62. if        0.86    x  lO-^-    (1500    -    12,100    x  0.07T 

x2  q0      0.035  0.07 

1    +    2.32    x   10 3  • 


8.OI4.  x  10 3 


By  trial  and  error,  we  can  solve  these  three  equations  with 

1    2       S 
three  unknowns,  x  ,  x  ,  and  —  .   The  results  are 

0 

qu 

x1  =  O.Olfl 
x2   =  0.051 

S  ft2 

=  0.05 . 


Calculating   the    cost   for   this    system,    we   obtain 


Y_      (APn)    =  3750    psi 


n=l 


3 


y    (APn)1/2  =  (iooo)1/2  +  (i25o)1/2  +  (i5oo)1/2 

n=l 


=   10 


5.7  /i/Psi 


4Pn  -  AP11"1        1000  250  250 

+  +  


,n-l 


0.035       O.Olfl       1.051 
=  3.957   x  10^  psi 
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IT   (Ren)2'8  =  14-.385  x  1012  . 
n=l 

Substituting  the  values  computed  into  the  cost  function, 

equation  (l\.2)  ,    we  have 

0.035  ,, 

Ct  =  O.7887  3.957  x  10^ 

0.035 
1  -  

0.07 
0.07  TO 

+  0.05  ( )    (0.4835  x  lcr1?  x  ij-385  x  1012 

0.07  -  0.035 

+  0.3266  x  10-7  x  3750  +  0.133  x  10"6  x  105.7) 
-  0.2286  x  10-8  x  1500 

=  3.4.192  x  io_5  —  =  0.284.8  ( : )  . 

lbm  1000  gal 

G.   Calculation  of  the  Cost  for  a   One-stage  System.   We 
have  the  total  number  of  independent  variables  of  this  system 
=2+1=3. 
(i)   We  choose  the  Reynolds  number  as  (34) 
Re  =  15,600 
the  operating  pressure  as 

A?   =  1000  psi 
and  the  discharged  concentration  as  (34-) 

x3  =  0.051  . 
Then  the  whole  system  will  be  fixed. 

Employing  equation  (5l)>  we  obtain 

1    "s   62.4 

0.051  1 (— ) 


0.035     q°   0.035 


100 


.     lo-^   (1000    -   12,100   x  0.05D 


=  1 


.  32    x   10 3 
(15,600) (/0 


.                   s       58.708I|.  s 
0.051  (28.57  -  (— )  =  1 

I  q°       1.025i4-      ' 


or 


1.0251].  s 

(28.57  -  19.607)  =  —  =  0.1565  . 

58.70  81+       q° 

The    cost   function   is 

0.035         1000 

Ct    =   O.7887    x   10_°   x  ( ) 

0.035      0.035 

1    -   

0.051 

0.051  /  17  p      D 

+  0.1565  x  ( )    (o.Lj.835  x  10"ir  x  (15,600)^-° 

0.051  -  0.035 

+   0.3266   x  lO"?    x   1000    +   0.133    x    (1000 )1/2    x   10"6) 

a           0.035 
-    0.2286   x  10"°  x   1000 

0.051    -    0.035 

=  3.9953  x  10~£  =  0.33261  ( )  . 

lbm  1000  gal 

(ii)  If  we  use  the  same  Reynolds  number  and  operating 
pressure,  but  use  the  discharging  concentration  at 

x1  =  0.07 
we  have 
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(      1  S      62.1+     0.86    x   10"^    (1000    -    12,100    x  0.07) 

0.07   '! }  =    1 

0.035       q°  0.035  0.07 

1   +    2.32    x   10 3  

1^.6563   x  10 3 


or 


/62.1+        0.86    x  10_i+  x   153 

(28.57  -  li+. 285)/  1 

I  0.035      1  +  0.031+9     J    q° 

or 

S    ll|..  285 

—  =  =  0.6302  . 

q°    22.668 

The  value  of  the  cost  function  is 

0.035    1000 
Ct  =  O.7887  x  10~b  x  ( ) 

0.035  0.035 

1  -  

0.07 

°'07         C  17  9  fl 

+  0.6302  ( )  (0.I+835  x  10"17  x  (I5,600)2*ti 

0.07  -  0.035   I 

+  0.3266  x  10~7  x  1000  +  0.133  x  (1000)1/2  x  10   j 

«    0.035 
-  0.2286  x  10"°  x  1000 

0.07  -  0.035 

=  6.355  x  10"^  =  0.529^  ( )   . 

lbm  1000  gal 

From  the  first  and  second  examples,  we  see  that  by  using  a 
gradually  increasing  operating  pressure  we  can  get  the  better 
result  (i.e.,  lower  cost)  than  by  using  the  same  operating  pres- 
sure throughout  the  system.   Prom  the  third  example  we  can 
recognize  that  our  proposed  sequential  process  is  better  than 
the  others.   Although  the  processes  presented  in  those  examples 
are  yet  to  be  optimized  completely,  we  believe  that  we  still  can 
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get  the  qualitative  conclusion  that  the  proposed  process 
better  than  other  processes  operated  under  these  similar  oper- 
ating conditions. 

Since  we  have  already  formulated  a  system  model  or  a  set 
of  system  equations  in  previous  chapters,  we  should  be  able  to 
use  the  model  to  find  the  optimum  design  and  operating  condi- 
tions of  the  system.   And  since  the  proposed  system  has  a  set 
of  well  defined  system  equations,  and  the  transformation  func- 
tions are  continuously  dif f erentiable  with  respect  to  the  state 
variables,  from  the  previous  two  parts,  we  see  that  the  system 
can  be  optimized  by  means  of  the  multi-level  approach  and/or 
the  discrete  maximum  principle.   The  main  difficulty  in  a 
numerical  solution  is  the  convergence  of  the  iteration  scheme. 
The  choice  of  the  step  size  factor,  k,  is  a  difficulty.   If  it 
is  too  large,  the  iteration  scheme  will  not  converge.   Yet,  if 
it  is  too  small,  convergence  is  extremely  slow.   The  step  size 
factor,  k,  can  only  be  adjusted  on  a  trial-and-error  basis  to 
achieve  convergence.   This  needs  a  lot  of  computer  time.   This 
work  will  be  left  for  future  work. 
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NOMENCLATURE  FOR  PART  THREE 

A,  An  =  Cross-section  area  normal  to  the  streamline  of  any 

membrane  separator  unit,  (ft2). 

K  .  d 

C      =  3.0£  x  10 •?  ■ constant.  ■ 

ScV3  Da 

CQ     =  Electrical  power  cost  ($/kw-hr). 
Cp    =  Pump  and  turbine  installed  cost  ($/kw). 
C      =  The  unit  cost  of  the  material  for  constructing  mem- 
brane separator  unit  ($/lbm) . 
C-j-     =  The  total  cost  per  pound  of  fresh  water  produced  ($/lbm). 
d     =  The  diameter  of  the  tubes  in  the  membrane  separator 

unit. 

sq  cm 

Da     =  Molecular  diffusion  coefficient  of  salt  ( ) . 

sec 
n 
E-,     =  Pump  work  of  high  pressure  pump  at  nth  stage. 

n 
Ep    =  Pump  work  of  circulation  pump  at  nth  stage. 

E?    =  Energy  recovery  from  the  blowdown  turbine  at  the  end 

of  the  process. 

f     =  Panning  friction  factor. 

Pn    =  The  volumetric  flow  rate  of  fresh  water  product  through 

ft3 

the  membrane  of  nth  stage  ( )  . 

ft2  -  hr 

pn    =   The  volumetric  flow  rate  of  salt  component  through  the 

ft-3 

membrane  of  nth  stage  ( )  . 

ft2  -  hr 

ft3 

K  ■   =  Membrane  constant  ( )  . 

ft2  -  hr  -  psi 
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m     =  The  total  number  of  tubes  within  each  membrane 

separator  unit. 
=  Total  number  of  stages  in  the  sequence  of  the  process. 
Pn    =  Pressure  within  membrane  separator  chamber  of  nth 

stage  (psi) . 
P     =  Atmosphere  pressure  ( li-j.-  7  psi). 
A?n   =  Pn  -  P  =  Pressure  difference  across  the  membrane  at 

nth  stage  (psi) . 
qn    =  Mass  flow  rate  of  brine  solution  discharged  from  nth 

stage  (lbm/hr) . 
q     =  Mass  flow  rate  of  feed  saline  water  (lbm/hr). 
R     =  The  radius  of  the  tubes  within  the  membrane  separator 

unit  =  1/2  d. 
Rn    =  Recycle  ratio  of  nth  stage. 
Ren   =  Average  Reynolds  number  at  nth  stage. 
S,  Sn  =  Membrane  area  of  one  membrane  separator  unit  (ft^). 
Sc     =  Schmidt  number. 

U     =  The  maximum  velocity  within  the  membrane  separator  unit 
u     =  Mean  velocity  within  the  membrane  separator  unit 

=  0.8  U  (for  turbulent  flow). 
un    =  Mean  velocity  within  the  membrane  separator  chamber 

of  nth  stage. 
v.„     =  Friction  velocity. 
Wn   .  =  Mass  flow  rate  of  fresh  water  produced  from  nth 

stage  (lbm/hr) . 
Wf    =  Total  mass  flow  rate  of  fresh  water  produced  from  the 

(ibm/hr)  system. 
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Wg     =  The  mass  of  the  she 11- and- tube  membrane  separator  unit 

of  nth  stage  (ibm). 
xn    =  The  stage  variable;  here  we  denote  the  mass  fraction  of 

the  salt  component  in  the  outlet  brine  solution  of  nth 

stage . 
xn    =  Average  mass  fraction  of  salt  concentration  within 

membrane  separator  chamber  of  nth  stage. 
Xq    =  Average  mass  fraction  of  salt  concentration  at  the 

membrane  surface  of  nth  stage, 
y     =  Distance  normal  to  membrane  boundary. 

Greek  Letters 

5  =  Momentum  transport  boundary-layer  thickness. 

6  =  Mass  transport  boundary-layer  thickness. 

\i      sq  cm 

=  Kinematic  viscosity  =  —  ( )  . 

f        sec 

\i  -   Viscosity  of  brine  solution. 

f  =  Density  of  brine  solution,  (lbm/ft3). 

fm  =   Density  of  material  of  constriction  (lbm/ft-^). 

0n  =  Decision  variable  of  nth  stage. 

7[n    =  The  osmotic  pressure  of  the  brine  solution  at  the 
s  ^ 

membrane  surface  of  nth  stage  (psi). 
Tq    =  Shearing  stress  at  the  wall. 
Ylf  =   Loss  factor. 

Y}m  =   Mechanical  efficiency. 
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>|0    =  Pump  efficiency, 
y)     =  Turbine  efficiency. 

=  Capitalization  charge  of  initial  cost  per  hour  in 
sure am  ($/hr) . 
c^    =  Allowable  stress  of  the  material  of  construction  (psi) 
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PROOF  FOR  THE  CONVERGENCE  OF  THE  PRICE- 
AD  JTJSTMENT  RULE 


Consider  the  nth.  subproblem.   For  some  prices  pn,  let  those 
components  of  Rn  which  are  zero  at  the  subproblem  extremum  be 
the  first  rn,  0  <  rn  1  /n.   That  is,  to  separate  Rn  >;  0  into 
equality  and  inequality  parts, 

R?On,  yn>N+1,  x""1)  =0,    i  =  1,  2,  ....  r"  (A. la) 

Rn(6n;  jA.H+1^  xn-l}  „  Q>        .  =  pn  rn  +  ^_  _ ;  ^  yn  _   (A_lb) 

Define  the  Lagrangian  for  the  nth  subproblem  as 

rn 
LTi   -  Sn  £  un  R1}  .  (A.1) 

j=l   J   J 

Substituting  the  expression  for  the  subobjective  function.  Sn, 

equation  (1.7),  Chapter  III,  Part  One,  into  this  equation,  the 

Lagrangian  becomes 

in  =   fn(en     yn}    +    (pn}T   T    (Qn^    xn-l}    _    (pn-l}T   xn-l 


+    (un)T   Rn(0n,    jn,    x11"1)     .  (A. 2) 

The  solutions  0n,  xn   ,  un  to  this  subproblem  satisfy  the 

following': 

<?Ln    9fn  m  2Tn       m  <^Rn 

-  +  (p^1  +  (unr  =  0  (A. 3) 

k  =  1,  2,  . . .,  wn 

m    ^Rn 


?Ln 

5fn                                  n,        c?Tn                                  m         ^R11' 

-J-    f    n ^  -1-                       (  ^\L 

n-1   ' 
^xk 

?x£  X                   ^x-1                   ^xg-1 

R1? •  =  0,    i  =  1,  2,  .  .  .,  r11   .  (A. 5) 
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3  following  independent  small  perturbation  of  pn  is 
made  at  each  subsystem, 

pn  _  pn  +  A pn  ^     n  =  1,  2,  . . . ,  N  ( A . 6 ) 

a  disturbance  will  alter  the  solutions  of  equations  (A. 3) 

.  (A. 5)  to 
3^-1  +Axn-1 

%*   +  AG11"1  (A. 7) 

and 

un  +  Aun  . 

The  variational  equations  can  be  obtained  from  equations 

(A. 3)  through  (A. 5)  as 

2fn                                       m      ^Tn                                       „      3Rn  n 

( )^    +    (p^+ApV    ( )^    +    (un+AuV    ( )A     =  0       (A.8)1 

k  k  k 

k  =   1 ,    2 ,     .  .  .  ,    w 

^xg-1^  °  ^xk 

^Rn 

+    (£n   +  AunL        ( -)&     =   0  (A.  9) 

*         ^xf1"1 
k 

xV       "~~        -1_   •  3  *     •     *     9  ^ 

(Rn)A   =  0  ,    i  =  1,  2,    . . . ,  rn  .  (A. 10) 

The  quantities  in  equations  (A.8)  through  (A. 10)  are  then 
expanded  in  the  Taylor  series  to  the  first-order  terms.   A 
representative  term  in  this  expansion  has  the  form 
?fn  Zfn        wn     22fn 


( — )A    =  —  +  H  &en, 

2  eg    <?eg   3=1  2e^aeg   J 


(   L  denotes  the  quantities  inside  the  bracket  and  is 
evaluated  at  the  perturbed  states. 


H5 


sn-l    ^2fn 

+   yr       


—  A  x1?"1  +  0(£2) 


(A. 11) 


J     k 
where  0(£^)  denotes  the  term  Including  second-order  and  those 

of  higher  order. 

When  all  terms  in  equation  (A. 8)  are  expanded  in  this 

manner,  we  obtain  the  following: 


•>  ,.i  i 


,.'  ,.ii 


in-  i. 


AQn.  +  y 


'k 


Ax. 

n-1  ftn  J 


n-_i_ 


dQji      j=i    3d1]  de^        J       j=i      dxLrx  e1,1 


s 


n. 


+    F      (E 
i=l 


n 


piJ   +  Apn 


/>!' 


dK 


w11  ^TJ.;i 

+  E      —  Ae" 


+  fc 


n-1 


j=i   JW^e. 

J        k 


Ax" 


j=i      ^j-^ek 


.  n 

+    II      (un  ■  +  Aun 
i=l 


n  n 

w 
i 


f  JR. 


k 


+  z 


d2Yd 

±—  ZiOn 

n.  p  Qn         j 


J=i    *e^.k 


,n-l 


3=1      ^x^-i^e" 


3= Ax 

n-1  ij  ftn  J 


n-1 


(A. 12) 


+  o(£2)    =  0    . 

A    similar   expression  is    obtained   by    expanding    equation    (A. 9). 
Expansion   of    equation    (A. 10)    gives 


,n 


>n 


R  •   +  Z      -AeVl 

3=1    ^©j         J       3=1      d*- 


n-1       ^nn 
s  ^Ri  n-1  ? 

^      A  x1?    1   +   0(£2)    =   0    .     (A. 13) 


>n. 


Since    each   of   the    first    ru  components    of   R.    is    zero    at    the 
original    price    Pq    ,    equation.    (A. 13)    reduces    to 
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E"  — -  ag1?  +  r   ^r  ax1:"1  +  o(£2)  =  o       U.iIO 

j=i  ^ex1   J   3=1   ^xn_1   J 
J     j  j 

Collecting  the  coefficients  of  A0.  and  x.    and  ignoring 

j  j 

all    terms    involving   the   product   of   two    or  more    increments, 
0(      ),    on   the    assumption  it   all    of   the   functions   fn(6    ,    J   ), 

Tn(6n,    xn_1),    and  Rn(6n,    yn,    xn_1)    depend    smoothly   on  x11"1   and 
0    ,    equation    (A. 12)    becomes 

-,,..n  n     /■     n2.nn 


Wn      (     ^? 

+  l-.    pv  — :  -'-  l.    u.-  — z  +  c. 


^e?      i=i"A£e?      i=i     1jef      j=i  ( ^en  <?e" 

k                             k                             k        J        v       j        k 
+    T      pn  —  +    Z      un  -(    A9n 

i=i    "^e^se^      i=i    1^ei?<^eg 

n-1  ,        ^2fn  sn  ,2* 

+     Z  n +     Z      P         TJ— 

*i      n    n-1  .,,  „n 


^       J  k  j 


k 


r-n  2  2Rn 


n  i        /     .     n-1 


+  T  un  — -±—  \  ax1: 

"-1  ,>Gr 
J 


1     IX'?"1   .'a'1    ! 


sn  „   ,?T?        rn  „  2R? 


APi~~  +  r  a^  — i-  =  o  .  (a. 15) 


i=i         -de?      i=i ~    -1  <?g£ 

k  k 


Prom   the    subob jective    function   S x,    equation    (2.5),    it    can  be 

seen   that 

^2Sn              02fn  „     ^2Tn 
_   ,     f    m  1 


^gx^g?     ze^eg  7     ^G'^eg 


K>  X1T7X  <a-i6> 
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32Sn  32fn  m         ^2Tn 


ax^^aef      Sx1?"1  ?ef  v    y     ^xn"13  9P 

J  k  j  k  j  k 

By  using    equations    (A. 3),  (A. 16),    and    (A. 17),    we    can    simplify 
equation    (A. 15)    to    obtain 


w£   r  *2sn        ,   xT  ^2Rn    )      _ 

i=i  i?enaen  *en  2en  J        J 

j        k  j        k 

n-1   A         ?2Gn                                    ,  2TDn       \ 
+    y      )  +   Anl1  I    £x 

3=1    Ux^^e?  ^-^e^  J 

^       J  k  j  k ' 

m     «?Tn  m     -9Rn 

+  -{<dpnY     +  &U11}     -  0    .  (A.  18; 

^  "  J      ion,   c   /   TQn 

An  analogous  development  applied  to  equation  (A. 9)  gives. 


w 


n. 


22Sn  „   ^2Rn 


T   — 7  +  <un)T  r   Zie1? 

j   K         j  ■  K 

sn-l     ^2sn  ^2Rn 


+  V 


"~3    "~k  ""j   '"k   ' 

+  {Apn)T  — —  +  4un>T  — -j  =  0  .              (A.  19) 

Now  we  define  the  following  vectors  and  matrices: 

<Qn>  =  <-^>  =  {q?}  -                                (a.  20) 

.X. 

-,   0        n.  n   -,        n  ,   n-1 

1=1,  2,  ...,  w  ,  w   +  1,  . . . ,  w   +  s 


22Sn       >2sn 

(Jn)  =  ( )  =  (— -)  (A. 21) 

lnn  _       _„n  o  n 

3  2R?        t  2R? 

:    =  ( — )  =  ( -)  (A. 22) 

k       k   J 


rn 


(Bn)  =  (Jn)  +  £"  V}   Kn  (A.  23) 


i=l   x  x 

to  n  ,  _n-l 

j  =  1,  2,     .  .  .  ,  w   +  s 

k  =  1,  2,  . . . ,  wn  +  sn_1 


n 


_L       _i_  «      5    •  •  *  9 

Note,  since  (Jn)  and  (K.)  are  symmetric,  the  matrix  (3n)  is 
symmetric  too.   Let 


and 


Ap"  -  0,    for    j  ±   i 


n-i 
Ap  .    =0,    for    all  j. 


Then  dividing  equations  (A.l8)  and  (A.  19)  by  Ap^1,  taking 

n 
tit  as  Ap.-  — >0,.and  using  the  definitions,  equations 

(A. 20)  through  (A. 23),  we  obtain 

^Qn       2Tn      ?Rn  m   „un 
(Bn)  {-U  =  -  <-±)  -  i  —  )1    {—}  (A.  21,) 


and 


3Rn   3Qn 


(— H)  <— )   =  0  (A. 25) 

n  =  1 ,  2 ,  .  .  .  ,  N 


Similarly,  letting 


.  n-1  t_   ~ 

APi    7=  0 


AP 


n-1 


and 


AP,  =  0, 
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-  0,    for    j  #  i 


for    all  j 


n-1 


dividing  equations  (A.l8)  and  (A.  19)  by  AP-r    anb  passing 


A  n-1 

Ap.  '  to  the  limit,  we  obtain. 


and 


(Bu) 


2Rn 

( ) 


#Qn 


^P 


n-1 


/Dwn+i 


<?Q 


n 


2Qn  Up?'1) 


=  0 


n  =  1,  2,  .  .  .,  N 


-1         -L  )  >     *  m     *     3  kJ 

3=1,  2,  ...,  rn 


n-1 


n 


( ) 

2Qn 


T 


<?u 


n 


^P?"1' 


,wn-fi 


(A. 26) 


(A. 27) 


is  the  (w   +  s    )  dimensional  unit  vector  with  a 

one  in  the  (wn+i)th  position. 

n  -1 
Premultiplying  equations  (A.2ij_)  and  (A.  26)  by  (B  )  .  with 

the  assumption  Bn  is  nonsingular,  and  substituting  the.  expres- 


,"n  /..  ^n 


>n 


.n-1 


sions    for  #  Q,  /c?  p  •       and     2  R  /dV^~         thus    obtained    into    equations 
(A. 25)    and    (A. 27),    we   have 


hRn  n       2Rn 

( )     (Bn)"       (■ 


T 


£u 


n 


2R 


n. 


n 


,n 


n. 


i  m 


2Q 


n 


)     (Bn) 


n,-l 


?T? 


2Q, 


n 


J 


(a. 28; 


J-  _Lj         *~  .?         *    •    *    9         ^ 
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and 


,:n  m ,     9un   i  3rh 


( )     .B")"1    ( )X( =    ( )     (Bn)_1   |Dwn+-       (     . 

i   =    1,    2,     ...,    s11"1 

r-    1 ,    2 ,     .  .  .  ,    N . 
For    simplicity,    we    define    the   following   expressions: 

(Cn)    =    ( )     (Bn)_1    ( y  (A. 30) 

<?Qn  3Q,n 

(CV)   =    ( )     (B11)-1!-^  (A. 3D 


^Pn|   =    (__)     (Bn}    1   £Dw   +ij  (A>32) 


2Rn 

re   Cn   is    an   rn  by  rn  matrix,-  G-?   is   an   rn-dimensional   column 

i 

n  n  n 

vector,  and  so  is  F. .   Note  that  since  B   is  symmetric,  C   is 

so  symmetric.   Thus  the  relations,  equations  (A. 28)  and 

. 29) ,  become 

(Cn)  J  -  -  (Cf  (A. 33) 


12        sn 


and 


°n>  (75^1  "  "  Ki  "■*> 


1 

i  =  1   2       sn-1 
Premultiplying  by  (Cn)~  ,  equations  (A. 33)  and  (A. 31+) 


become 


*u* 

-   (cn)' 

1V\  '  " 

^U11       s 

1  - 

=  (cn) 

^pflj 
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1  [c^l  (A.  35) 

and 

{^|         '  (A.  36) 

n.  -1 
Premultlpiying  by  (B  )    and  making  use  of  equations 

(A. 35)  and  (A. 36),  equations  (A.2I4.)  and  (A. 26)  become 

•  ^Qn^  1  ,3T^         1      £Rn  m       -,  ,  „  , 

{— r  =  -  (B")"1  f  —if  +  (B^r1  ( )   (Cn)_1  /(£/      (A. 37) 

-L  5         J  •   •   •   a      O 

] =  (Bn)  1  [Dw  +1   -  (Bn)  X  ( )L    (Cn)    W     (A. 38) 

Now,  in.  order  to  obtain  the  elements  of  —  ,    we  return  to 
the  definition  of  the  vector  of  excess  demand  defined  as  (see 
equation  (I.76,  Chapter  III,  Part  One), 

En(P)  =   xn  -  Tn(6n,  x11"1). 

Taking  the  first  partial  derivatives  of  En(P)  with  respect  to 

p  ,  yie±ds 

?En    5xn    ?Tn(Sn,  x""1) 

=  — (A. 39) 

n.  =  1,  2,     .  .  .  ,    N 
m  =  1 ,  2,     .  .  . ,  N . 


Recall  that  P  =.  (p  .  p  :  ...  :  p  )  is  a  matrix  and  pn, 
n  =  1,  2,     .  .  .  ,    N,  Is  a.  column  vector. 


•essary  conditions  for*  the  subproblem  to  be  e 

tremized,  equations  (A. 3),  (A. 4),  and  (A. $)  ,    it  can  be  recog- 

Lzed  th    bhe  nth.  subproblem  solutions  xn   ,  %n ,    and  un  are  the 

notions  of  pn  and  pn   .   And  since  the  solutions  xn  are  de- 

Lned  by  the  (n+l)th  subsystem,  and  they  are  only  dependent 

on  the  values  of  pn  and  prj   ,  and  not  dependent  on  the  value  of 

0En 
Thus  the  only  nonvanishing  terms  of  are  those  where 

m  =  n-1,    n,  and  n+1.   They  are: 

,?sn    ?xn   STnan-'\  Hn)  ?Tnan"\  en) 


n  =  '2, ;  3,  ■...,  N 

*En    3S?1   fT^S?1"1,  8n) 


(A. ij.0) 


and 


?  pn      ^  p"n 

n  =  1,  2,  .  .  .  ,  N 


fl-rpn-1    „^n-l    ornn-l,.^n-2  ^n-l->    ^^n-1 
^E       dx  2T         (x    ,  6    )    2x 


(A. 1+1) 


3  Pn       ?n  ^Pn  #  Pn 


(A. I(.2) 


n  =  2,  3,  ...,  N. 

By  using  the  definition  of  Qn,  equation  (A. 20),  the  ele- 
ments of  equations  (A.40),  (A. 41),  and  (A. 42)  become 

9*1  ,«£,*  ,  5Qn 


(  9Qn)        l   ,?nn  1; 


3-9.  i&q:xj    '  <?p.- 

k  =  1,  2,  .  . 


,n 


i  =    2       sn_1 
n  =  2,  3,  .  .  .  ,  N, 


123 


and 


^n  ,  o^,n+l 


9vl 


[ 


2T 


^P? 


?Qn 


k 


1,    2, 


2E 


n-1 
k 


<?P 


n 


i  =   1,    2, 


n.  =   1,    2, 


=      |Dw"+kjT     ^ 


,n 


,n 


,    N   -    1, 


<2P? 


.n-1 


(*Q 


n 


P 


n 


.n 


(a. i)4; 


(A. 1+5) 


k  =   1,    2,     ... 

1  J-  ^         *—  y  •    *     •    y  " 

n  =2,    3,    ..'.,    N. 

Substituting   equations    (A. 37)    and    (A. 38)    into    equations 
A.ljJ,     (A.i|lj_)>    an^    (A.ljJ?),    we    obtain 


2E 


^P 


n-1 


^k)T 


^Q 


n 


r1  {Dwll+i} 


^ 


,m-l 


2Rn   T 


'   — -f         (Bu)    "    ( ) 

!■£    —     A. ,      c. }      ,  .  .  f      S 

i   =   1.    2 s11"1 


n,-l 


Cu) 


Ki 


U.l;6) 


2E° 

k 

2P? 


=      D 


n 


k  =   2,    3,     .  .  .  ,    N 


n+1 


w  +kl   t    /^n+lx-1    iVw        +i 


(B^1)^    j 


L> 


,n+l 


Dwn+k|T    (Bn+irl    {^_Z^    (cn+l}-l 


^Q 


n+l 


{  — ^  /        (Bn)     X 


<?Q 


n 


PT 


n 


9Ql 


n 


pn+l 
i 


12J+ 


-      —  (B^)"1    ( )T    (Cn)_1   (g?  (A. 47) 

l^on  J  sop  l    ij 


and 


k  =  1,  2,  .  .  .,  sn 
i  =  1,  2,  ...,  sn 
n  =  1,    2,    ...,    N   -    1, 

i 


/D»nV    (B")"1^ 


1  )  l*Qn 


<?T 


fDvjn+k|T    (Bn}-l    ( )T    (cn)-l/Gn)  u>[j_8) 

k  =  1,    2,    .  .  . ,    sn_1 


X 


=•12  sn 


n  =   2,    3,    .  .  .,    N. 

n        n  n    . 

3y    substituting   the  matrices   C    ,    G^,    and  F-^    into    equations 

(A.lj.6),     (A.ij.7),    and    (A.i;8),    and  noting   that   Cn   and  Bn   are    sym- 
metric,   it   can  be    seen   that 

B^l  3E^  *e£       3^1       JE^1       JET1 


n-1        „„n-l    '    „    n 


,~n  i  «n 


> 


JV±  ^Pk  ^Pi       9Vk         3Vll  ipg 

and 

^En      ^E11"1  m 
( -)  =  ( )X  . 

dvn~L  ^pn 

?En     ^E 

That  means  that  and  —  are  symmetric  matrices.   This  enables 

d  Pm     3  V 

dE 
us  to  consider  the  entire  matrix  of  —  rather  than  just  its  sym- 

metric  part. 


1^5 


BE 
To  prove  that  is  negative  definite,  let  us  introduce  a 

quadratic  form  defined  as 


N 
n=l 

where 


,6  =  t   {Hnj  T  (B^)"1  (Hn|  (A.Lj.9) 


Hn  =  Bn  ( )  {?vn}  +  Bn  {A""1  (A. 50) 

^Pn  <?pn_1       J 

n  =  1,  2,     .  .  .  ,    N, 

and  i  Xn  j  is  a  column  vector  with  elements  of  sn  arbitrary  real 
numbers.   When  equations  (A.I4.9)  and  (A.^O)  are  combined,  the 
quadratic  form  becomes 

N      f  <?Qn  ^  Qn  T 

p  -  r  6n  (— )  {^1 +  Bn  ( — )  h11"1})  (Bn)_i 

n=l   l   X  5pn      C  ^P11"1  ■     7 

•     (Bn    (  — )    {xnj    +   Bn    (— r-T)-  (An_1 


^Pn  #p 


n-1 


£   ffx        (~ >    (B  }    (— ;: 

n==l     I  dp  d  P 


+    fx11"1)         ( -)       (Bn)       ( )    {xnj 

+  £x*)T  (— )T  (bV  (-^-)  (a-1) 

^pn  ?pn-!      C  j 


/V-1}  T  (-^-)T  (b-)t  (^-)  fx^/r    .  (a. 51: 

1  J         •  ^P11"1  2  V11-1 


Since  (B)  is  a  symmetric  matrix  and 


a  scalar  quantity,  their  transpose  should  be  identical, 
is 

(Bn)T  =    (Bn)     , 


Lt      >  oPn  ,?pn-l      <•  J  J 


=  An-ll*    (Z^L)*    (Bn)T    ffi    fx*/      .  (A. 52) 

<■  J         ^pn-1  apn  J 


Thus    equation    (A. 51)    becomes 


T     (M      (  —  )T   (Bn)T    (— )    (x* 
n=l    *«  ^pn  7Pn 


+   2{xn-l}*    (Z*L,*    (Bn}T    (ftf|     £x„) 


,n  ^n 


+   Un_1         ( r)T    (bV    ( -)    (A""1         .  (A. 


53) 


( r)T    (bV    (  — 

Is    is   a    quadratic   form   of  X's,    which  can   be   written   as 

p   =  (X}T    (A)    |X]  (A. 54) 


where 

) 
l  s    s 


fx}T=({xl}T(x2jT...     {X»1TJ 


N 

y  sn 

n=l 

dimensional  row  vector. 

We  shall  now  prove  that 
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(a)  =  —    .  (a. 55: 

Equation  (A. 55)  will  be  proved  if  we  can  show  that  the  co- 

T  0EU 

efficient  of  {/Vn)   f^j  i-n   P  ^s  an&   t^-*3  coefficient  of 


\XU}T   l^"1]  in  p  is 


^pn-l 


(a)   Proof  that  the  coefficient  of  (\n)      fxn  1|  in  (3  is 

^En 
equal  to  .   Prom  equation  (A.5l),  "the  coefficient  of 

(Xn]T  [An_1}  term  is 

2Qn  j>Qn 

( )   (.B11)1  ( 


3Vn  ^Pn-1 

Taking  the  transpose  of  the  equation  (A.  21).)  we  have 

Bn)  ( )    =  ( )X  (F/V 

^pn  i      <?pn 

-  -  ( )  -  ( )  ( )  .  (A. 56) 

Postmultiplying  equation  (A.  5&)  "by 

( T} 

we  obtain 

( )T  (B11)1  ( -) 

<?pn       2pn_1 


=    -    ( )  ( -)  -  ( )   ( )  ( -)   .        (A. 57) 

lon  (A. 2?)  is  substituted  into  equation  (A. 57),  we 
obtain 

( r  (bV  ; — )  =  -  ( — )  ( -)     .  i.$Q) 

This  equation  is  exactly  tbe  equation  needed  to  complete 

z'm   proof. 

i  T  • 
(b)   Proof  that  the  coefficient  of  jxn/   j  Xnl  in  (3  is  equal 

5En  J 


to 


Prom  equation  (A.5l),  the  coefficient  of  fxnj   ( \n(    term 

is 

( )T  (Bn)T  ( )   and   ( r  (bV  ( )   . 

2Pn         5Pn         ^Pn  <?Pn 

By  postmultiplying  eouation  (A. 56)  by  ( ),  we  obtain 

( )X  (bV  ( ) 

d Pn  2  Pn 

J1,Tn   ^Qn     aun    ?Rn   ^Qn 

=  -  ( )  ( )  -  ( )   ( )  ( )   .         (A. 59) 

<?Qn  dvTi  3vn        <?Qn   ^Pn 

By  using  the  condition,  equation  (A. 25),  equation  (A. 59)  becomes 

( )   (Bn)   ( )  =   -    ( )  ( )  .  (A. 60) 

<?Pn  ^pn       ^Qn   ^pn 

Taking  the  transpose  of  the  equation  (A. 26),  we  obtain 
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^Qn+1  v  m     m   ?un+1  m  5Rn+1 


( r  (bV  =  (D)T  -  ( )X  ( -)  (A. 61) 

where  (D)  is  a  (w    +  s  )  by  s.  matrix,  with  as   by  s   unit 

matrix  I  under  a  wn   by  s"n  null  matrix. 

5Qn+1 

Postmultiplying  by  ( ),  equation  (A. 61  becomes 

<?pn 

( r  (bV  ( ) 


o 


^n  o  v^n 


=    (D)T    ( )    -    ( )T    ( -)     ( )       .  (A. 62) 

2pn  <?Pn  3Qln+1  <?Vn 

By  using   the    condition  of    equation    (A. 27),    equation    (A. 62) 

becomes 

2Qn+1   T  T     ^Qn  T      ^Qn" 


'Bn)~    ( )    =    (D)~    ( )       .  (A. 63) 


3  pn  ^  pn  ^  pn 

Combining   equations  (A.60)    and    (A. 63),    the    coefficient    of 
{\n}X   {Xn)term   is 

(D)T    ( )  -    ( )     ( )       .  (A.6^) 

^Pn  ^Qn       JVn 

Comparison  with  equation  (A.ijij.)  shows  that  this  quantity  is 
precisely  . 


9Vn 

Thus  equation  (A. 55)  is  proved.   This  means  that  —  is  the 

<?p 
ma  trice  of  the  quadratic  form  (3.   But  from  the  definition  of  (3, 

equation  (A.lj.9),  the  quadratic  form  is  a  sum  of  smaller  quad- 

ratic  forms,  \Enr       (Bn)    {_En]  .      Thus  we  link  the  negative 

def initeness  of — —  with  the  negative  definiteness  of  the  ma- 

trices  (Bn)   . 


»n\  /cn>  _1 


nee  (B  ) (B  )   '  =  (I),  which  is  positive  definite,  the 
gative       Lteness  of  the  matrices  (Bn)  will  guarantee  the 
;ative  definiteness  of  the  ma  trice,  (Bn) 

By  the  definition,  equation  (A. 23),  the  matrices  B   are 
related  to  the  matrices  of  second  partial  derivatives  of  the 
subobjective  functions  S   and  of  the  constraints  R. .   The 
structure  of  the  matrices  B  will,  of  course,  vary  with  P,  since 

Lch  constraints  are  active  and  which  are  not  depends  upon  P. 
This  leads  to  the  following  theorem  (6) . 

3 or em.   If  for  all  n,  the  subobjective  function  Sn  and 
all  the  constraints  Rn  are  concave  (for  maximization  problems) 
in  the  arguments  x    and  Gn  for  all  real  values  of  ?,  and  if 
at  least  one  of  these  functions  is  strictly  concave,  the  price- 

justment  rule,  equation  (1.17),  Chapter  III,  Part  One,  is 
asymptotically  stable  in  the  large,  and  convergence  to  P  is 

monotone  in  |E   . 

n      n 
Proof.   3y  the  hypotheses,  the  matrices  J  and  K.  m  equa- 
tions (A. 21)  and  (A. 22)  are  negative  semidef inite,  with  at  least 
one  negative  definite,  for  ail  P,  0  ,  and  xJ   ,  n  =  1,  2,     . . . , 
N   :  =  n   2       vn 

J-i  %  -*-        -j-  *       9  *  •  *  9  * 

Then  by  equation  (A. 23)  all  Bn,  and  hence  (Pn)   ,  are 
jative  semidef inite,  with  at  least  one  negative  definite. 
.is  implies  that  in  the  expression,  equation  (A.Lj_9),  the  quad- 
ratic form  [3  is  negative  definite  for  all  \    and  P,  which  guar- 

2E 
antees  the  negative  definite  of  the  matrice  —  . 
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It  Is  shown  in  section.  2,    Chapter  III,  Part  One,  that 

asymptotic  stability  in.  the  large  of  price-ad justment  rule, 

dV 

equation  (1.17),  Chapter  III,  Part  One,  requires  that  —  be 

dt 

negative  definite  for  all  P  or,  eouivalently,  that  in  the 

expression  (1.20),  Chapter  III, Part  One,  —  should  be  negative 
definite  for  all  P.   Thus  the  proof  is  completed. 
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A    comparative  and  critical  study  of  the  multi-level  system 
theory  and  the  maximum  principle  was  carried  out.   While  the 
two-level  structure  of  the  multi-level  theory  was  proved  to  be 
identical  to  the  discrete  maximum  principle  for  simple  recycle 
processes,  it- appears  to  be  extremely  difficult,  if  not  impos- 
sible, to  equate  the  multi-level  theory  to  the  discrete  maximum 
principle  for> systems  which  are  more  complex  than  the  two-level 
structure.   However,  it  is  plausible  that  we  can  develop  an 
optimization  technique  in  which  both  the  multi-level  theory  and 
the  discrete  maximum  principle  can  be  jointly  used.   In.  an 
attempt  to  develop  such  a  method  the  maximum  principle  was 
extended  to  systems  with  inequality  constraints  by  using  the 
Kuhn  and  Tucker  complementary  slackness  principle  which  is  one 
of  the  tools  employed  in  developing  the  multi-level  theory. 
The  system  model  and  performance  equations  of  the  reverse- 
osmosis  water  -purification  process  were  developed  for  the  pur- 
pose of  optimizing  the  process  by  means  of  the  multi-level 
approach  and/or  the  discrete  maximum  principle. 


